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Abstract 

This  thesis  is  a  continuation  of  a  previous  effort  which  developed  a  finite  element 
solution  of  Schrodinger's  Equation.  The  purpose  of  this  research  is  to  extend  this 
previous  work,  and  develop  a  chemical  laser  engineering  tool  for  the  identification  of 
transition  lines.  Identification  of  laser  transitions  for  a  new  chemical  gain  medium 
requires  knowledge  of  Einstein’s  Coefficients.  These  transitions  rates  can  be  obtained  by 
solving  Schrodinger's  Equation  for  diatomic  molecules  using  the  method  of  finite 
elements.  Experimental  vibrational  eigenvalues  for  a  given  electronic  state  are  used  to 
determine  the  molecular  potential  surface  which  yields  the  closest  numerical  result.  A 
non-linear  minimization  routine  is  used  to  hunt  for  this  surface  by  adjusting  parameters  of 
energy  functions  such  as  the  Harmonic,  Morse,  Lennard-Jones,  and  Mie  potentials.  For 
each  set  of  new  parameters  selected  by  the  minimization  routine,  the  method  of  finite 
elements  is  used  to  solve  Schrodinger's  Equation  .  The  eigenvalues  from  these  solutions 
are  then  compared  to  the  experimental  values.  Through  this  iterative  process,  the  best 
potential  surface  is  isolated.  Franck-Condon  factors,  which  are  proportional  to  Einstein’s 
coefficients,  can  be  computed  with  the  numerical  eigenfunctions  from  two  different 


potential  surfaces  found  in  this  way. 


This  numerical  technique  was  able  to  isolate  potential  surfaces  whose  eigenvalue 
solutions  had  relative  errors  better  than  10'3  and  10'6  percent  when  compared  to  the 
analytical  solutions  of  the  Harmonic  and  Morse  oscillators,  respectively.  Comparisons  of 
the  wavefunctions  also  yielded  excellent  agreement.  Initial  work  with  H2  (X  'Sg+) 
verifies  the  lower  eigenstates  can  be  approximated  by  the  Morse  potential  with  an 
anharmonicity  term  of  1 .0912  inverse  a.u.  and  a  dissociation  energy  of  0.177  Hartrees. 


IDENTIFICATION  OF  MOLECULAR  LASER  TRANSITIONS  USING  THE 


FINITE  ELEMENT  METHOD 


I  .Introduction 


Background 

The  Air  Force  has  on-going  research  to  develop  new  and/or  better  battlefield  laser 
weapon  systems  which  will  prove  powerful  and  economical  in  the  field.  This  thesis 
research  project  contributes  to  this  effort  by  aiding  the  identification  of  molecular 
candidates  for  chemical  laser  systems  in  a  cheap  and  effective  manner.  Military 
applications  of  these  laser  systems  include  countermeasures  against  infrared  guided 
missiles,  theater  defense,  weapons  guidance,  and  other  battlefield  applications. 

To  determine  whether  a  molecule  is  a  good  candidate  for  lasing,  the  molecule’s 
energy  transition  rates  must  be  understood.  These  transition  rates,  calculated  using 
Einstein’s  Coefficients,  identify  if  a  sufficient  population  inversion  can  be  achieved  to 
establish  lasing  between  two  energy  levels  of  an  atom  or  molecule  (1:179-183,616-624). 
Often,  this  data  is  not  available  or  incomplete  leaving  the  researcher  unable  to  analyze  the 
molecule  for  its  lasing  potential.  At  this  point,  the  investigator  must  endeavor  in  a  time 
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consuming  and  possibly  expensive  laboratory  research  to  obtain  this  data,  and  then 
evaluate  the  molecule  as  a  lasing  candidate. 

However,  Franck-Condon  probability  factors  which  are  related  to  these  transition 
rates  can  be  obtained  through  theoretical  calculations  if  sufficient  spectroscopic  data  is 
available.  Commonly,  molecular  spectroscopists  use  the  semi-classical  Rydberg-Klein- 
Rees  (RKR)  method  with  an  Inverted  Perturbation  Approach  (IP A)  refinement  for  these 
calculations.  The  Air  Force  Institute  of  Technology  in  the  early  1980’s  (2;  3;  4)  used  a 
RKR-IPA  model  developed  by  Vidal  and  Scheingraber  (5)  as  a  research  tool  for  previous 
work  related  to  this  thesis. 

The  RKR  method  uses  constants  experimentally  derived  from  spectroscopic  data 
to  calculate  the  classical  turning  points  of  vibrating  diatomic  molecules.  These  turning 
points  are  then  used  to  construct  the  potential  energy  curve  and,  thereby,  gamer  the 
wavefunctions  needed  to  calculate  the  transition  probabilities.  RKR  does  not  rely  on  the 
actual  spectral  data  nor  the  energy  eigenvalues  for  this  operation,  but  molecular  constants. 

Problems  can  arise  from  reducing  large  sets  of  spectroscopic  data  into  molecular 
constants.  “Because  of  the  inadequacy  of  the  Bom-Oppenheimer  separation  of  the  total 
molecular  energy  into  electronic,  vibrational,  and  rotational  parts,  a  large  number  of 
molecular  constants  must  be  introduced  to  account  for  the  energy  level  structure  of  the 
molecule.  These  constants  appear  in  the  expressions  for  the  levels  in  an  often  complex 
and  nonlinear  manner  so  that  their  determination  poses  a  burdensome  problem  ... 
Moreover,  additional  difficulties  may  arise  from  the  need  to  determine  the  molecular 
constants  using  proper  statistics.”  (6:38)  With  these  difficulties  in  mind,  and  possible 
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systematic  errors  in  the  rotational  constants  “within  the  experimental  limit  (are)  sufficient 
to  cause  substantial  non-physical  behavior  in  the  potential  obtained  on  inversion.”  (7:248) 

This  non-physical  behavior  often  manifests  itself  in  the  repulsive  branch  of  the 
potential.  The  repulsive  arm  can  either  have  a  ripple  or  may  bend  over  (7:244)  making 
the  potential  double-valued  and  non-physical,  typically  at  the  higher  energy  levels  which 
have  greater  experimental  uncertainty  or  are  unknown.  This  situation  leaves  the 
transition  calculations  for  the  higher  energy  levels  suspect.  Tellinghuisen  and  Henderson 
(8)  circumvent  this  problem  by  replacing  the  repulsive  branch  of  the  potential  with  the 
Morse  potential  (9).  However,  this  technique  does  not  attack  nor  change  the  fundamental 
problem  of  with  the  RKR-IPA  method,  that  being  this  method  sometimes  yields 
pathological  potential  energy  curves. 

In  1984,  Shankland,  Dorko,  and  Ostdiek  developed  a  new  approach  to  find  the 
potential  energy  curve  for  diatomic  molecules  (10).  Their  method,  unlike  RKR,  is  a 
quantum  mechanical  approach,  and  uses  the  actual  absorption  and  emission  spectral  data 
for  its  calculations.  This  numerical  technique,  involves  a  finite  element  solution  of 
Schrodinger's  Equation  for  a  given  potential  energy  curve.  Each  computed  eigenvalue, 
along  with  its  associated  experimental  eigenvalue,  is  then  compared  and  a  residual  is 
recorded.  The  potential  energy  function  is  constructed  so  that  a  non-linear  minimization 
routine  can  then  vary  the  parameters  within  this  function.  Then  for  each  new  potential 
curve  selected  by  the  minimization  routine,  Schrodinger's  Equation  is  re-solved.  Thereby, 
through  an  iterative  process,  the  potential  curve  which  yields  the  smallest  residual  is 
found.  The  wave  functions  can  now  be  garnered  from  this  potential.  Along  with  the 


1-3 


wave  functions  from  another  electronic  state,  the  transition  probabilities  are  then 
calculated. 

The  Shankland-Dorko-Ostdiek  (SDO)  method  has  several  advantages  over  the 
RKR  method.  One,  it  is  completely  a  quantum  mechanical  approach,  not  a  hybrid  of 
classical  and  quantum  mechanics.  Therefore,  its  solutions  are  rooted  in  physical  law  and 
theory  as  understood  today.  Second,  it  bypasses  the  problems  RKR  encounters  with 
molecular  constants  by  using  eigenvalues  derived  directly  from  the  spectroscopic  data. 
This  spectroscopic  data  does  not  have  to  be  reduced  into  spectroscopic  constants,  but  only 
assigned  to  the  proper  transition.  Lastly,  any  potential  function  can  be  used  for  the  fit. 
Possibilities  include  the  Morse,  Lennard-Jones,  Mie  potentials,  or  even  a  custom  potential 
function.  Therefore,  this  method  will  never  render  a  potential  function  which  is 
pathological.  The  accumulative  effect  of  these  advantages  means  that  this  method  has  the 
prospect  of  yielding  results  of  greater  accuracy  than  RKR,  especially  at  the  higher  energy 
states. 

However,  the  SDO  method,  after  initial  development,  underwent  very  little  testing 
and  evaluation.  Many  questions  went  left  unanswered.  Will  this  numerical  technique 
converge  for  sophisticated  potentials  such  as  a  Morse  or  a  Mie?  Will  this  method  work 
for  a  real  molecule? 
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Objectives 


The  main  objective  of  this  research  is  to  validate  the  SDO  method  as  a  reliable 
numerical  solution  of  Schrodinger's  Equation.  Reliability  and  accuracy  will  be 
determined  by  comparing  the  numerical  eigenvalues  and  functions  with  the  analytical 
solutions  of  oscillators  such  as  the  Harmonic  and  Morse.  Second,  demonstrate  that  this 
method  can  be  used  for  finding  the  potential  energy  curve  of  a  molecule  which  is 
previously  known. 

Scope 

The  scope  of  this  research  project  is  to  validate  the  use  of  this  numerical  technique 
to  predict  energy  transitions  and  their  associated  probabilities  for  diatomic  molecules. 
Validation  will  include  comparisons  to  the  analytical  solutions  of  a  simple  harmonic, 
analytical  solutions  of  an  anharmonic  oscillator,  and  one  real  molecule  with  a  large, 
accurate,  empirical  knowledge  base.  This  study  is  only  concerned  with  the  analysis  of 
diatomic  molecules,  even  though  this  code  with  some  modification  could  be  extended  to 
linear  triatomics  or  even  more  complicated  molecules.  No  attempt  will  be  made  to 
collect  any  spectroscopic  data  in  a  laboratory  environment,  only  published  spectroscopic 
data  will  be  used.  Also,  no  effort  will  be  made  to  apply  this  code  to  predict  undiscovered 
data  for  some  arbitrary  molecule. 
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Approach 


This  thesis  report  is  organized  in  a  linear  fashion  starting  with  theory  and 
background,  and  ending  with  validation  studies.  The  objective  of  this  organization  is 
convince  the  reader  that  this  numerical  approach  is  a  reliable  tool  for  calculating 
transition  probabilities  for  diatomic  molecules. 

Chapter  II  begins  to  lay  the  theoretical  foundation  by  relating  Einstein’s 
Coefficients  to  the  Franck-Condon  Factors  through  a  quantum  mechanical  treatment. 
This  development  is  followed  by  a  discussion  of  molecular  forces,  and  the  electronic 
shepherding  which  is  this  binding  force.  These  molecular  forces  can  be  approximated  by 
the  Harmonic  and  Morse  Oscillators.  The  analytical  solutions  of  these  potentials  are 
shown  in  prelude  to  a  discussion  about  the  Franck-Condon  Principle.  Unfortunately,  real 
molecules  do  not  behave  like  these  potentials,  which  leads  to  the  requirement  for 
numerical  techniques  such  as  the  SDO  method.  The  rest  of  the  chapter  is  dedicated  to 
this  method  of  finite  elements  as  used  to  solve  Schrodinger's  Equation. 

Chapter  III  continues  with  this  thought  by  explaining  how  the  SDO  method  is 
implemented  on  a  computer.  Discussions  include  how  spectroscopic  data  is  converted 
into  a  format  acceptable  to  the  model,  how  the  model  searches  and  finds  the  “best” 
potential  surface  with  a  non-linear  minimization  routine,  and  how  the  Franck-Condon 
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Factors  are  calculated  using  the  eigenfunctions  garnered  from  this  optimal  surface.  These 
discussions  include  a  basic  recipe  for  executing  each  section  of  the  code. 

Finally,  chapter  IV  demonstrates  that  this  numerical  technique  is  a  valid  approach 
for  solving  Schrodinger's  Equation.  This  validation  is  accomplished  by  comparing  the 
analytical  and  numerical  solutions  of  the  Harmonic  and  Morse  Oscillators.  Then,  a 
known  potential  surface  for  the  hydrogen  molecule  is  used  as  a  test  to  see  if  indeed  this 
potential  surface  could  find  this  potential  in  the  blind. 

A  concluding  chapter  summarizes  these  results  and  presents  recommendations  for 
more  research. 
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II .  Literature  Background  and  Theory 


Introduction 

The  motivating  force  behind  this  research  is  to  develop  a  tool  a  laser  engineer  can 
use  to  analyze  the  electronic  and  vibrational  transitions  taking  place  within  the  gain 
medium.  In  order  to  analyze  these  transitions,  the  quantum  mechanical  or  molecular 
properties  of  this  gain  medium  must  be  known. 

Therefore,  this  chapter  introduces  the  theory  and  concepts  needed  to  understand 
how  the  Shankland-Dorko-Ostdiek  (SDO)  numerical  technique  reveals  these  quantum 
mechanical  properties  of  the  molecule.  First,  an  overview  of  transition  theory  lays  the 
foundation  for  what  knowledge  is  required  of  the  molecule  for  the  laser  engineering 
analysis  to  happen,  namely  the  Franck-Condon  Factors  related  through  Einstein’s 
Coefficients.  To  calculate  these  factors,  quantum  mechanics  is  invoked  to  gamer  the 
molecule’s  wave  functions.  Therefore,  a  short  discussion  follows  which  sets  up 
Schrodinger's  Equation,  a  solution  of  which  gives  the  wave  functions,  and  shows  how 
several  different  potential  functions  can  be  used  to  approximate  the  molecular  forces. 

The  SDO  technique,  a  solution  for  Schrodinger's  Equation  numerically  using  a  finite 
element  method,  is  briefly  explained.  Finally,  with  the  wave  functions  in  hand,  the  last 
section  outlines  the  procedure  for  using  the  wavefunctions  to  calculate  the  Franck- 
Condon  Factors,  and  describes  the  principle  on  which  they  are  founded. 


Laser  Transitions 


A  laser  (Light  Amplification  by  Stimulated  Emission  of  Radiation)  is  a  device 
which  utilizes  the  natural  energy  transitions  of  atoms  and  molecules  to  produce  a 
coherent,  amplified,  and  monochromatic  source  of  electromagnetic  radiation.  For  an  in- 
depth  discussion  of  lasers  see  Verdeyen  (1). 

Briefly,  however,  a  fundamental  requirement  to  initiate  lasing  is  to  establish  a 
population  inversion  between  two  transition  levels.  This  non-equilibrium  state  is 
accomplished  by  pumping  the  molecules  of  the  gain  medium  from  a  lower  energy  state  to 
an  upper  state  with  an  external  energy  source.  These  excited  molecules  then  relax  back  to 
their  ground  state,  giving  up  their  energy  through  either  radiative  or  non-radiative  events. 
The  wavelength  of  the  radiation  is  determined  by  the  difference  in  energy  between  the 
two  states  on  which  the  transition  happened.  The  rate  of  the  radiative  processes  described 
by  the  Einstein  Coefficients,  coupled  with  the  non-radiative  rates  determine  if  the 
population  inversion  can  be  sustained.  If  the  population  inversion  is  lost,  the  lasing 
activity  ceases.  Figure  II- 1  illustrates  the  transition  process  for  a  hypothetical  four  level 
laser. 
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Figure  II-l  This  schematic  diagram  shows  a  possible  arrangement  of  the  energy  levels  for  the  gain 
medium  of  a  laser.  The  pumping  action  promotes  some  of  the  ground  state’s  population  (labeled  0) 
to  the  highest  energy  level.  From  this  energy  level,  the  excited  population  relaxes  back  to  the  ground 
via  the  path  and  rates  as  indicated.  Note:  the  rates  are  also  proportional  to  the  population  of  the 
states  initiating  the  transition.  Lasing  actually  occurs  on  the  transition  from  state  2  to  1,  where  the 
field  stimulates  emission  adding  constructively  more  energy  to  the  field.  Einstein’s  Coefficients  are 
labeled  A  and  B  (for  a  more  detailed  definition  see  the  following  discussion),  the  rate  of  non-radiative 
processes  are  indicated  by  the  symbol  K,  I  is  the  intensity  of  the  lasing  field,  and  c  is  the  speed  of 
light. 


The  determination  of  these  process  rates  called  Einstein’s  Coefficients,  therefore, 
are  critical  to  the  transition  analysis  of  any  laser.  These  coefficients  can  be  derived 
through  time-dependent  perturbation  theory  of  Schrodinger's  Equation  (2:TPl-TP12b; 
3:267-277;  6:12-18,24-25).  If  the  Hamiltonian,  a  quantum  mechanical  representation  of 
the  atom’s  or  molecule’s  energy,  is  perturbed  by  an  in-coming  time-dependent 
electromagnetic  field  then  the  rate  of  transition,  or  absorption  for  this  example,  is  found 
to  be 
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where  R  is  the  rate  of  transition  from  state  k  to  state  m,  k  is  the  initial  state  of  the  system, 
m  is  the  final  state,  c  is  the  speed  of  light,  h  is  Planck’s  constant*,  s0  is  the  permittivity  of 
free  space,  I  is  the  intensity  of  the  electromagnetic  field  at  the  frequency  comk 


corresponding  to  the  transition,  p  is  the  electric  dipole  operator,  and 


2 


is 


defined  as  the  transition  moment. 

In  Figure  II- 1,  for  example,  k  would  be  state  1  and  m  would  be  state  2.  Therefore, 
the  absorption  rate  from  energy  state  1  to  2,  indicated  by  Bi2  I/c,  must  be  the  same 
quantity  as  defined  by  Equation  II- 1 .  With  that,  setting  them  equal  to  each  other  and 
keeping  the  units  straight  yields 
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Now  if  the  pump  was  turned  off,  and  the  system  allowed  to  reach  thermal 
equilibrium,  the  relationship  governing  the  populations  of  states  2  and  1  must  still  hold  at 
thermodynamic  equilibrium.  Each  radiative  process  going  down  in  energy  must  equal 


*  Special  note:  due  to  a  font  limitation,  hbar  will  always  be  written  as  h/2n.  Also,  the  2n  factor  may  be 
simplified  into  the  expression. 
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each  process  going  up.  Therefore,  the  time  rate  of  change  of  each  population  for  level  1 
and  2  must  be  equal  to  each  other.  Mathematically  this  means 
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must  be  true.  Where  Ni^  is  the  population  of  each  respective  level,  p(cc>2i)  is  the  energy 
density  of  the  field  previously  referred  to  as  I(comk)/c,  A21  is  the  probability  per  second  for 
spontaneous  emission  of  a  photon,  B12  p((02i)  is  the  probability  per  second  for  absorption 
of  a  photon,  and  B21  p(o)2i)  is  the  probability  per  second  of  spontaneously  emitting  a 
photon.  Since,  the  system  is  at  thermal  equilibrium,  p(©2i)  must  be  equal  to  energy 
density  defined  by  Planck’s  Blackbody  formula.  This  formula  is 
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where  V21  is  the  frequency  of  light  equal  to  <x>2\/2n,  k  is  Boltzman’s  constant,  and  T  is  the 
temperature  at  thermal  equilibrium  (6:  24). 


Solving  Equation  II-3  for  p(©2i)  and  setting  it  equal  to  Planck’s  Blackbody 
formula,  “Einstein  forced  the  fit  with  identification  of  various  interrelationships  between 
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the  coefficients”  (4: 1 82).  Designating  the  degeneracy  of  each  respective  level  as  gis2 , 
Einstein  found 
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Casting  in  general  terms,  Einstein’s  Coefficients  become 
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The  transition  moment  can  be  further  simplified  for  the  simple  case  of  diatomic 
molecules.  Using  the  Bom-Oppenheimer  (5  )  approximation,  the  wavefunction  for  each 
state,  m  and  k,  can  be  separated  into  independent  electronic,  vibrational,  and  rotational 
functions.  Therefore,  the  total  diatomic  wavefunction  for  each  state  can  be  expressed  as 
(6:136-137) 
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where  T'ei  represents  the  electronic  wave  function,  Xv  represents  the  vibrational  wave 
function  of  the  molecule,  Xj  represents  the  rotational  wavefunction  of  the  molecule,  and  r 
and  R  represent  the  sets  of  electronic  and  nuclear  coordinates,  respectively. 

Since  changes  in  electronic,  vibrational,  and  rotational  states  are  possible  both  r 
and  R  dependencies  must  be  considered  in  the  dipole  operator  p  when  calculating  electric 
dipole  (El)  transition  probabilities.  Therefore,  the  dipole  operator  may  be  written  as 


E  =  -  Z  er-  +  E  eZNR  N  +  Vnucl 
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Where  e  is  the  charge  of  an  electron,  i  is  an  index  for  each  electron  in  the  molecule,  Z  is 
the  charge  of  the  nucleus,  and  N  is  the  index  for  each  nuclei. 

Ignoring  the  rotational  behavior  of  the  molecule  (for  simplicity)  the  transition 
moment  becomes 
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Each  dipole  operator  only  operates  on  its  corresponding  wavefunction.  Therefore,  each 
dipole  operator  is  distributed,  and  after  some  factoring,  each  then  operate  on  its 
corresponding  wavefunction.  The  following  simplifications  are  then  made.  Immediately, 
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the  second  term  goes  to  zero  because  the  two  electronic  wavefunctions  are  orthogonal  to 
each  other.  The  first  term,  however,  does  not  go  to  zero  because  the  two  vibrational 
wavefunctions  are  from  two  distinct  basis  sets  and  not  necessarily  orthogonal.  The 

integration  of  (V™  \xel  'f'Jj  'j  over  r  yields  the  transition  moment,  Me  (R) .  Me  (R)  is  a 

function  of  nuclear  position  R  because  the  electronic  wave  functions  depend  on  both  r 
and  R.  However,  if  Me  (R)  is  assumed  to  vary  slowly  with  R,  then  further  simplification 
yields 


(m| (x | k)  =  Me(R)  (x:|x||) 

=  Me(R)  (v'|v") 


(II- 10) 


where  Me(R)  is  the  average  value  of  the  electronic  transition  moment.  From  here  on,  a 
shorthand  is  introduced  to  represent  the  upper  electronic  state’s  vibrational  wavefunction 
as  v' ,  and  v"  will  represent  the  lower  electronic  state’s  vibrational  wavefunction.  The 
probability  of  the  electronic  transition  is  then  proportional  to  the  total  transition  moment 

squared,  or  in  mathematical  terms  (|m  p  .  Then,  substituting  Equation  II- 10  into 
Equation  II-6  produces  the  final  result 
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where 


is  defined  as  the  Franck-Condon  Factor.  The  Franck-Condon  Principle 


is  discussed  in  detail  later.  Thus,  Einstein’s  Coefficients  are  directly  proportional  to  the 
Franck-Condon  Factors  for  diatomic  molecules. 


Diatomic  Molecules 

Clearly,  the  Franck-Condon  Factors  are  important  to  laser  engineering  in  regards 
to  transition  analysis.  To  calculate  these  factors,  Schrodinger's  Equation  must  be  solved 
for  the  molecule  of  interest.  The  following  discussion  outlines  the  theory  and  procedure 
for  solving  Schrodinger's  Equation  analytically. 

For  simplicity,  imagine  a  reduced  mass  representation  of  a  vibrating  molecule  that 
is  fixed  in  space,  i.e.  not  translating,  and  not  rotating.  Then  Schrodinger's  Equation  may 
be  expressed  as  (7:98-101) 
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where  H  is  the  Hamiltonian  Operator,  E  is  the  energy  of  molecule,  'Fvib  is  the 
wavefunction  which  describes  the  nuclear  position,  p  is  the  reduced  mass  of  the 
molecule,  R  is  the  intemuclear  separation,  and  V(R)  is  the  potential  energy  of  the  system 
as  a  function  of  R.  If  the  potential,  V(R) ,  is  known,  then  Equation  11-12  written  in  matrix 
form  can  be  diagonalized,  thereby,  gamering  the  wavefunctions  and  energies  and 
ultimately  the  Franck-Condon  Factors  are  calculated  from  the  wavefunctions. 

Physically,  V,R)  represent  fte  force,  related  by  FA„B  -  ,  which  binds 

the  two  atoms  together  into  one  molecule.  This  force  arises  from  the  variation  of  total 
electronic  energy  with  intemuclear  separation.  In  other  words,  the  molecular  electron 
cloud  shepherds  the  two  nuclei  keeping  them  bundled  together.  This  shepherding  process 
overcomes  the  nuclear-nuclear  repulsive  forces  yearning  to  dissociate  the  molecule. 

This  shepherding  process  is  extremely  complex.  In  principle,  however,  the 
complete  Hamiltonian  could  be  written  down  for  the  entire  molecular  system.  This 
expression  would  include  all  of  the  forces,  such  as  the  Coulombic  electron-nuclear 
attractive  and  electron-electron  repulsive  forces,  the  rotational  terms,  spin  coupling, 
relativistic  corrections,  and  so  on.  But  these  Hamiltonians  are  extremely  difficult  to 
solve,  even  for  the  simplest  of  molecules  like  H2.  Therefore,  one  way  to  avoid  this 
problem  is  to  lump  all  of  these  terms  together,  call  it  V(R) ,  and,  thereby,  analyze  only  the 
accumulative  molecular  potential. 

This  potential  energy  surface,  however,  cannot  be  represented  by  any  one  general 
function  for  all  diatomic  molecules.  There  are,  however,  some  general  properties  the 
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potential  must  possess  to  model  the  physical  behavior  correctly.  One,  it  must  be  smooth 
and  continuous.  This  condition  ensures  the  force  is  always  finite.  Second,  the  function 
needs  to  be  always  single  valued.  That  is,  at  any  given  intemuclear  separation,  only  one 
energy  configuration  is  allowed.  Next,  as  intemuclear  distance  goes  to  infinity  the 
molecular  forces  must  go  to  zero  to  allow  for  dissociation,  therefore,  the  slope  of  the 
potential  must  go  to  zero.  Conversely,  as  the  intemuclear  distance  approaches  zero  the 
nuclear-nuclear  repulsive  forces  increase  dramatically.  Therefore,  the  slope  of  the 
potential  needs  to  be  steep  and  negative  to  model  this  force.  Finally,  a  stable  molecule 
implies  an  equilibrium  position  exists  where  the  force  is  zero,  i.e.  the  potential  energy  is 
at  a  minimum.  Figure  II-2  illustrates  such  potential  energy  surface. 
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Intermolecular  Force  Model 


V(R) 


Figure  II-2.  This  figure  illustrates  the  intermolecular  force  model.  The  potential  energy  curve,  V(R), 
drawn  with  the  solid  line  exemplifies  all  of  the  characteristics  needed  to  model  the  intermolecular 
forces,  FAoB,  between  atoms  A  and  B.  The  steep  repulsive  branch  captures  the  large  forces  the  two 
nuclei  experience  at  small  R,  and  the  attractive  branch  on  the  right  allows  the  molecule  to  dissociate 
at  large  internuclear  separations.  The  well  defines  the  region  of  a  lower  energy  configuration 
centered  by  the  radius  of  equilibrium,  Re,  where  the  energy  is  a  minimum.  De  is  dissociation  energy 
as  defined  from  the  bottom  of  the  well  to  the  dissociation  limit.  The  graph  of  the  harmonic  oscillator 
shows  where  this  approximation  is  valid. 


A  Taylor  expansion  of  some  general  potential,  V(R) ,  which  meets  these 
conditions,  about  the  equilibrium  position,  Re,  yields 
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where  the  value  of  the  potential  at  R«  is  -De  (the  dissociation  energy).  The  first  derivative 
vanishes  at  Re  because  the  potential  is  at  a  minimum.  A  constant,  k,  is  the  curvature  of 
the  potential.  Figure  II-2  illustrates  this  Taylor  expansion  as  the  dashed  curve.  This 
truncated  expansion  is  immediately  recognized  as  the  harmonic  oscillator. 

The  solution  of  Schrodinger's  Equation  for  the  harmonic  oscillator  is  (8:176-202) 


(11-14) 


Where  En  is  the  energy  of  state  n,  w  is  the  angular  frequency  of  the  oscillation  equal  to 


—  ,  k  is  the  force  spring  constant,  p  is  the  reduced  mass,  'Tn  is  the  eigenfunction  of 
In¬ 


state  n,  a  -  oi p  —  ,  q  =  va  x ,  x  =  (R-Re),  and  Hn(q)  is  the  Hermite  polynomial 
^  h 


Figure  II-3  shows  4  wave  functions  for  this  simple  oscillator. 
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n  =  0  n  =  1 


Figure  II-3.  This  graphic  shows  four  eigenstates  of  the  harmonic  oscillator  (9:67-73).  Each 
eigenfunction  is  plotted  at  its  corresponding  energy  level  for  illustrative  purposes  only,  the  amplitude 
of  the  waves  have  units  of  probability/(unit  length)I/2.  Notice  each  wavefunction  is  symmetric  about 
the  radius  of  equilibrium,  Re,  and  the  number  of  nodes  corresponds  to  its  state  number.  The 
amplitude  squared  of  these  wavefunctions  is  the  probability  of  finding  the  reduced  mass  particle  at 
that  location.  Notice,  the  wavefunctions  imply  the  particle  can  actually  tunnel  through  its  classical 
turning  point,  where  the  wavefunctions  crosses  the  potential  curve.  Plotted  with  Mathematica®. 
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The  utility  of  this  harmonic  approximation  is  limited  to  the  region  very  near 
equilibrium.  Morse  proposed  instead  the  following  function  as  a  suitable  approximation 
for  the  diatomic  molecular  energy  curve  (10) 


V(r)  =  D  e  2a(R  Re)  -  2Dee  a(R  Re) 


(11-15) 


where  a  is  referred  to  as  the  anharmonicity  term.  This  variable  has  units  of  inverse 
length.  See  Figure  II-4  for  an  illustration  of  the  Morse  potential. 

The  width  of  the  Morse  potential  at  half  max,  AW1/2,  is  inversely  proportional  to 
a.  Therefore,  the  larger  a  is,  the  stronger  the  molecular  forces  are,  and  the  narrower  the 
well  becomes.  The  dissociation  energy  De  stretches  the  well  vertically.  De  and  a 
determine  the  number  of  bound  states  which  are  allowed  by  this  potential.  Large  De  and 
small  a  yield  large  numbers  of  bound  states. 
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V(R) 

Morse  Anharmonic  Oscillator 


Figure  II-4.  This  figure  illustrates  the  Morse  Potential.  This  curve  exhibits  all  of  the  characteristics 
required  to  model  molecular  forces.  In  fact,  this  model  is  a  good  approximation,  however,  no 
molecule  is  known  to  exactly  mimic  the  Morse  Oscillator.  The  width  at  half  max,  AW1/2,  is  inversely 
proportional  to  a.  Therefore,  the  larger  a  is,  the  stronger  the  molecular  forces  are,  and  the  narrower 
the  well  becomes.  The  dissociation  energy  parameter  De  stretches  the  well  vertically. 


Morse  analytically  calculated  the  eigenvalues  and  functions  for  this  oscillator  as 
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where  n  is  a  non-negative  integer  which  labels  the  eigenstate,  to  =  a 
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Associated  Laguerre  Polynomial  (1 1 :  59-60;  12:  725-726)  given  by 
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Morse  commented  that  this  was  the  first  solution  of  Schrodinger's  Equation  which 
yielded  a  finite  number  of  discrete  energy  states.  In  fact,  the  condition  on  b  (see  Equation 
11-17)  for  the  Laguerre  Polynomial  to  be  defined  implies  n  must  always  be  less  than  d. 
Expressed  mathematically  this  condition  is 
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n  =  0  n  -  3 


Figure  II-5.  This  figure  illustrates  4  eigenfunctions  of  a  Morse  Oscillator.  Each  eigenfunction  is 
plotted  at  its  corresponding  energy  level  for  illustrative  purposes  only,  the  amplitude  of  the  waves 
have  units  of  probability/(unit  length)172.  Notice  each  wavefunction  is  asymmetric  about  the  radius  of 
equilibrium,  R*.  R*  is  denoted  by  the  dashed  vertical  line.  The  wavefunctions  are  stretched  towards 
the  classical  turning  point  associated  with  the  attractive  branch.  The  large  amplitudes  imply  the 
probability  of  finding  the  reduced  mass  particle  at  this  location  is  greatest  here.  Notice  the  ground 
state  (n=0)  is  very  similar  to  the  Harmonic  Oscillator’s  (see  Figure  II-3).  This  region  of  the  well  can 
be  approximated  by  a  harmonic  oscillator.  Plotted  with  Mathematica®. 
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Circling  back,  the  Morse  Oscillator  can  itself  be  approximated  by  the  Harmonic 
Oscillator  near  equilibrium.  The  force  constant  k  of  the  Harmonic  Oscillator  can  be 
related  to  the  Morse  Oscillator  by  expanding  it  into  a  Taylor  Series(see  Equation  11-13). 
This  relationship  is  found  by  setting  the  second  order  term  of  the  expansion  equal  to  the 
Harmonic  Oscillator  yielding 

k  =  2De  a2  (11-19) 


Franck-Condon  Principle 

“The  Franck-Condon  Principle  governs  the  intensity  of  spectral  transitions 
between  the  vibrational  levels  of  different  electronic  states  of  molecules.”(13:78-81). 
Therefore,  as  discussed  earlier,  the  Franck-Condon  Principle  is  central  to  the  process  of 
determining  laser  line  transition  probabilities.  The  following  section  describes  the 
Franck-Condon  Principle  in  more  detail. 

By  recognizing  the  great  difference  between  the  masses  of  the  electrons  and 
nuclei,  the  Bom-Oppenheimer  approximation  can  be  invoked.  Consider  two  different 
molecular  states,  then  each  state  can  be  separated  into  electronic,  vibrational,  and 
rotational  wavefunctions.  Recalling  Equation  II-7,  the  wavefunction  were  separated  and 
written  as 
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m  =  T.: 


Then  Equation  II-9  allows  us  to  calculate  the  El  transition  between  these  states  as 
(ignoring  rotation) 
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k)=  Mc(i?)(X^|Xk 


=  Me(R)  v' 


where  Me  (R)  is  the  average  electronic  transition  moment,  and  { v' 


(11-21) 


v")  is  defined  as  the 


vibrational  overlap  integral.  The  double  prime  notation  traditionally  represents  the  lower 
electronic  state,  and  single  prime  the  upper  state.  The  electronic  transition  moment,  a 
resultant  of  the  two  electronic  wavefunctions  and  the  dipole  operator,  is  averaged  over  the 
range  of  R  and  is  considered  constant.  The  probability  of  making  a  transition  from  v’  to 
v”  is  then  given  by  the  absolute  value  of  Equation  11-21  squared  where  the  Franck- 
Condon  Factor  is  defined  as  the  quantity 


qv'v"  — 
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,  where  0<  qvV,  <  1 

(11-22) 


where  v’  and  v”  are  two  vibrational  wavefunctions  each  from  a  different  electronic  state. 
For  illustrative  purposes,  Figure  II-6  shows  several  Morse  wavefunctions  depicting 
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Franck-Condon  Principle 


Figure  II-6.  This  figure  illustrates  the  concept  of  the  Franck-Condon  principle.  Transitions  occur 
vertically  where  the  overlap  integral  is  non-zero.  Following  the  Born-Oppenheimer  approximation, 
the  electronic  transition  happens  so  fast  that  the  nuclei  can  be  considered  to  be  at  rest  In  other 
words,  from  the  start  of  the  transition  to  the  end  of  the  transition,  the  nuclei  will  not  have  moved. 
Therefore,  a  transition  can  only  occur  from  one  state  to  the  next  if  the  nuclei  have  a  probability  of 
existing  at  that  location  in  both  states. 


the  vibrational  states  v’  and  v”.  These  vibronic  (a  contraction  of  the  words  vibrational 
and  electronic)  transitions  seldom  occur  without  a  change  in  the  vibrational  state  label. 
For  example,  the  upper  electronic  state  might  be  at  vibrational  level  3  and  then  transition 
down  to  the  lower  electronic  state  with  a  vibrational  state  of  1  as  depicted  in  Figure  II-6. 

These  transitions  are  considered  to  be  vertical.  That  is,  the  transitions  take  place 
at  a  fixed  nuclear  separation.  The  Bom-Oppenheimer  approximation,  based  on  the  large 
differences  between  the  masses  of  the  electrons  and  nuclei,  establishes  also  that  the 
relative  kinetic  energy  of  the  electrons  is  extremely  large  compared  to  the  nuclei. 
Therefore,  during  a  vibronic  transition  the  electronic  molecular  cloud  reconfigures  itself 
so  quickly  that  the  nuclei  virtually  haven’t  moved  in  this  time  frame.  With  the  nuclei 
being  fixed  in  position,  the  probability  of  this  vibronic  transition  can  only  be  non-zero  if 
and  only  if  the  nuclei  have  a  probability  of  existing  at  this  location  in  both  states. 

Summing  all  of  the  Franck-Condon  Factors  from  one  vibrational  state  in  the  upper 
electronic  state  to  all  of  the  vibrational  states  of  the  lower  electronic  state 
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where  closure  has  been  invoked  to  eliminate  the  summation.  Since  v’  is  orthogonal  with 
itself,  the  expression  reduces  to  1 .  Or  simply  put,  a  single  vibrational  state  in  the  upper 
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electronic  level  has  a  100%  probability  of  transitioning  to  some  vibrational  level  in  the 


lower  electronic  state. 

Certainly  then,  if  the  vibrational  wavefimctions  are  known  for  each  electronic 
state,  each  and  every  Franck-Condon  Factor  can  be  calculated. 


Numerical  Solution  of  Schrodinger's  Equation 

Numerical  methods  are  essential  for  any  realistic  molecular  potential  functions. 
The  Harmonic  and  Morse  potentials  are  functions  which  can  be  reduced  analytically. 
Unfortunately,  nature  does  not  provide  molecular  potential  curves  which  behave  so  nicely 
as  these.  Real  molecular  potential  functions  can’t  be  modeled  by  simple  functions  due  to 
the  sophistication  and  complexity  of  the  electronic  molecular  cloud.  This  complexity 
means  an  analytic  solution  for  the  vibronic  wavefunction  may  never  be  found  for 
molecules. 

Previous  to  1984,  the  Air  Force  Institute  of  Technology  used  the  customary  tool, 
the  RKR-IPA  method,  to  numerically  calculate  the  Franck-Condon  Factors  of  diatomic 
molecules.  Even  with  fine  tuning  and  “work  arounds”,  the  RKR  method  was  limited  to 
the  lowest  vibrational  states  (See  Chapter  I,  The  Introduction).  Therefore,  Shankland, 
Dorko,  and  Ostdiek  sought  to  develop  a  new  and  better  method.  The  technique  they 
developed,  which  exploits  the  method  of  finite  elements,  abandoned  all  of  the  methods 
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previously  employed  by  the  spectroscopy  community.  Even  though  this  method  of  finite 
elements  was  new  to  this  community,  it  was  certainly  not  to  the  engineering  world. 

The  emergence  of  computers  was  the  genesis  of  the  finite  element  method. 
Computers  enabled  scientists  and  engineers  to  tackle  problems  previously  unsolvable. 
Mathematical  techniques  such  as  the  classical  Rayleigh-Ritz  method,  variational  calculus, 
and  Galerkin’s  weighted  residuals  method,  much  developed  at  the  turn  of  the  century, 
were  assembled  into  what’s  now  called  the  finite  element  method.  This  method  was  first 
developed  in  the  early  1950’s  “to  solve  continuum  problems  in  elasticity  using  small 
discrete  elements  to  describe  the  overall  behavior  of  simple  elastic  bars...”  (14:1-2)  Over 
time,  this  method  grew  in  sophistication  and  has  been  used  successfully  to  analyze  an 
eclectic  assortment  of  problems,  such  as  electromagnetism,  heat  conduction,  fluid  flow, 
and  mechanical  stresses,  even  in  three  dimensions. 

The  SDO  method  finds  the  “best”  empirical  molecular  potential  in  an  exhaustive 
search.  Central  to  this  search  algorithm  is  a  numerical  solution  of  Schrodinger's  Equation 
solved  with  the  method  of  finite  elements.  The  following  discussion  explains  briefly  how 
the  finite  elements  method  is  used  to  solve  Schrodinger's  Equation.  Later,  the  Computer 
Modeling  Chapter  explains  how  this  technique  is  used  in  an  iterative  fashion  to  find  the 
best  potential. 

Schrodinger's  Equation  in  the  coordinate  representation  is 

(H-E)Tvil  =  --L-V2T  U  V(R)-EIJ'FIlb  =  0  (11-24) 

8ti2p  1 
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where  V(R)  caputures  only  the  vibrational  energies.  For  simplicity  the  rotational  energies 
are  ignored.  This  simplification  means  that  this  algorithm  will  not  handle  rotational 
energies.  Rewriting  Equation  11-24  in  the  time  independent  Dirac  notation  yields 
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(H-25) 


The  wavefunction  is  assumed  to  be  spherically  symmetric,  therefore,  in  the 


P2  h2  1  d  n2  d 

coordinate  representation - > - —  R  — 

2  p  8tt2  p  R2  <3R  dR 


and  integration  is  conducted 


over  the  volume  element  ck =4tiR23R  .  With  the  substitution  of  vFvjb  =  #  Ostdiek 

R 

showed  (15:18-36),  with  the  aid  of  integration  of  parts,  that  Equation  11-25  can  be  cast  in 
the  weak  form  in  the  coordinate  representation  as 
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where  now  dr  =  3R .  U  is  now  treated  as  the  spherical  wavefunction. 

To  employ  the  finite  element  method,  a  grid  is  laid  over  the  space  of  integration. 
The  goal  of  finite  elements  is  to  approximate  the  integral  (the  value  of  the  function)  in 
each  discretized  element.  Since  this  problem  has  been  reduced  to  one  dimension,  each 
element  has  the  step  length  of  h  (not  to  be  confused  with  Planck’s  constant).  Each  local 
element  is  then  transformed  by  a  natural  coordinate  system  to  ease  integration  later. 
Figure  II-7  illustrates  the  finite  element  grid. 


Finite  Element  Grid 


Natural  Coordinate  System 


hlj  - 

h - *■ 

hl2 

R  i-l  R  ;  R  i-2 


Figure  II-7.  This  diagram  depicts  the  Finite  Element  Grid  and  Natural  Coordinate  System.  The 
wavefunction  and  potential’s  value  and  slope  are  approximated  at  each  gridpoint  R{.  The  stepsize  of 
each  element  is  h.  The  stepsize  does  not  necessarily  have  to  be  constant,  but  for  the  algorithm 
developed  here,  h  is  constant.  The  natural  coordinates  Ij  and  12  are  introduced  to  ease  the  integration 
later  on. 
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The  following  relationships  are  true  for  the  natural  coordinate  system 


R  —  Rj  +  l2h, 

R  =  Rj  j  -  l,h, 

1+1  1  (11-27) 

h  =  Ri+1  -  Rj , 

1.  +12  =1 

where  Rj  is  the  location  of  the  ith  grid  point,  and  1)  and  h  are  the  natural  coordinates. 

The  finite  element  method  converges  as  long  as  the  requirements  for  completeness 
and  compatibility  are  satisfied  (16:1 14-1 15;15:22).  Namely,  the  value  of  U  and  its  slope 
must  be  continuous  at  each  grid  point.  Therefore,  U  can  be  approximated  at  the  local 
element  to  be 


Ue(R)  =  Ue(l,,l2)  =  U0  f,(l,,l2)  +  U0  f2 (1. » 12)  +  Uj  f3(l,,l2)  +  U,  f4(l„l2)  (11-28) 

where  the  superscript  e  means  this  formula  is  only  valid  for  the  local  element  in  question. 
Uo  and  Ui  are  constants  which  have  values  of  wavefunction  located  at  the  endpoints  Rj 

r  t 

andRj+i.  U0  and  U,  are  constants  which  represent  the  slope  of  U  at  Rj  and  Rj+i, 

respectively.  The  interpolating  polynomials,  fj-4,  turn  off  and  on  in  order  for  Equation  II- 
28  and  the  derivative  of  this  equation  to  meet  boundary  conditions.  For  example, 

Ue(R=Rj)  should  equal  Uo,  and  Ue  |r  =  Ri+1j  should  equal  Uj  ,  etc. 
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Since,  V(R)  is  required  for  the  integration  of  Equation  11-26,  in  a  similar  fashion 
as  this,  it  is  also  approximated. 

The  natural  coordinates  are  then  used  to  build  the  interpolating  functions  as  cubic 
polynomials.  The  basis  set  for  cubics  in  1|  and  b  is  {h  ;  li  b ;  li  b  ;  b  }•  The  cubic 
interpolating  functions  are  then  substituted  into  Equation  11-26.  Now  the  integrals  are 

t 

simple  polynomials  completely  in  terms  of  b  and  b,  and  the  constants  Uo,  U0  ,  Ui,  and 

f 

U,  .  Using  the  formula  (17:37;  15:25) 
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computing  the  integrals  becomes  straight  forward  and  trivial.  As  an  example,  the  first 
integral  in  the  numerator  of  Equation  11-26  (after  some  clever  manipulation)  evaluates  to 
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(11-30) 


In  a  similar  manner,  the  rest  of  Equation  11-26  is  tackled.  The  second  integral  in  the 
numerator,  which  includes  the  potential  energy,  breaks  apart  into  four  separate  integrals. 
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Therefore,  including  the  denominator,  six  different  integrals  are  evaluated  and 
manipulated  into  forms  similar  to  Equation  11-30.  For  complete  details  see  reference  15. 
An  important  note,  these  matrices  are  a  consequence  of  the  finite  element  method,  and 
should  not  be  confused  with  matrices  found  in  the  formalization  of  quantum  mechanics. 

Now  that  each  of  the  six  integrals  have  been  transformed  into  matrices  for  the 
local  element  only,  it’s  now  time  to  assemble  the  global  system.  This  equation  is  the 
matrix  system  which  represents  the  entire  grid,  not  just  the  local  element.  Each  local 
element  has  a  neighboring  element  on  either  side.  Therefore,  the  boundaiy  values  of  each 
element  (its  edges)  should  match  its  neighbor.  With  this  in  mind,  a  global  matrix  can  be 
assembled  by  overlapping  the  matrices  on  top  of  each  other.  See  Figure  II-8  to  see  how 
this  is  done. 
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Figure  II-8.  This  figure  illustrates  how  the  global  finite  element  matrix  is  assembled.  Each  grid 
element  consists  of  two  endpoints.  At  each  endpoint,  the  unknowns  are  the  wavefunction  value  and 
its  slope.  Hence,  that’s  four  unknowns.  But,  each  element  has  a  neighbor  to  each  side,  and  at  these 
common  boundaries  their  functional  values  need  to  match.  Therefore,  each  local  matrix  overlaps 
with  its  neighbors  matrix.  See  Equation  11-30  for  a  sample  matrix  equation  which  could  represent 
the  element.  This  diagram,  however,  only  depicts  the  first  two  elements  being  overlapped.  The 
building  process  actually  continues  until  every  grid  element  has  been  included.  Notice  the  square 
matrix  is  banded,  it  consists  of  one  main  diagonal  plus  three  lower  and  upper  diagonals.  Therefore, 
the  vast  majority  of  its  elements  are  zero.  The  “X”  in  the  diagram  indicates  matrix  multiplication. 
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For  simplicity,  Figure  II-8  only  shows  the  global  matrix  being  built  for  one 
integral.  Remember,  there  are  five  other  integrals  to  do.  Five  integrals  are  summed 
together  in  the  numerator  and  manipulated  into  a  form  similar  to  Equation  11-30  which  is 
just  one  of  these  five  integrals.  The  denominator  is  fashioned  in  a  similar  way. 

If  we  identify  the  column  vector  (this  contains  the  unknown  U  and  U’  values)  as 
Z,  the  square  matrix  in  the  numerator  as  H,  and  the  square  matrix  in  the  denominator  as  S 
Equation  11-26  becomes 


(11-31) 


In  summary,  Z  contains  information  about  the  wavefunctions,  H  describes  the  energy, 
and  S  contains  information  about  the  grid.  The  S  matrix  is  a  normalization  operator. 

To  find  the  eigenvalues  and  eigenvectors,  the  H  matrix  needs  to  be  diagonalized. 
However,  this  form  does  not  lend  itself  immediately  to  an  eigenvalue  problem. 
Therefore,  Equation  II-3 1  is  manipulated  into 
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HZ 


ESZ  =  0 
ELLtZ  =  0 


HZ 

h(lt)~'ltz  -  elltz  =  o 

(11-32) 

L",h(lt)''ltZ  -  EL'‘LLtZ  =  0 
XY  -  EY  =  0 

where  X  =  L_1H  |lt  j  and  Y  =  LT  Z ,  and  LLT  is  the  Cholesky  decomposition  of  S 
(15:32-35). 

X  has  the  same  eigenvalues  as  H.  Therefore,  the  eigenvalues  can  be  immediately 
found  by  diagonalizing  the  matrix  X.  The  diagonalization  is  easily  done  with  “off  the 
shelf’  math  routines  (see  Chapter  3).  Recovery  of  the  eigenvectors,  on  the  other  hand, 
requires  an  extra  step.  The  eigenvectors  are  simply  related  to  Y,  that  is,  Z  =  (l/f’Y.  The 
eigenvectors  are  normalized  by  the  relationship  N  =  Z  SZ,  where  N  is  the  normalization 
factor. 

The  rank  of  the  matrices  goes  as  m  =  2  (ne+ 1),  where  m  is  the  rank  of  the  matrix, 
and  ne  is  the  number  of  elements.  Therefore,  when  X  is  diagonalized  m  eigenvalues  will 
be  found.  Also,  Z  (which  contains  the  eigenvectors)  is  an  mxm  matrix.  Each  column  of 
Z  is  an  eigenvector,  where  each  column  contains  (ne+1)  pairs  of  data.  This  data  is  the 
eigenfunction’s  value  and  slope  at  each  gridpoint. 


11-32 


The  accuracy  of  this  solution  is  related  to  the  number  of  elements.  As  the  number 
of  elements  increases  the  stepsize  decreases  (assuming  the  endpoints  of  the  grid  are 
fixed).  The  convergence  is  (14:37-38) 


e(x) 


<  ch 


n+1 


(11-33) 


where 


e(x)  is  the  residual  between  the  exact  solution  and  the  numerical  solution,  c  is  a 


constant,  h  is  the  stepsize,  and  n  is  the  order  of  the  interpolation  polynomial.  Hence,  for 
cubic  elements,  the  error  decreases  as  order  h4. 

An  interesting  consequence  of  this  solution,  is  that  the  number  of  elements  must 
be  increased  to  decrease  the  stepsize  h  in  order  to  increase  the  accuracy  of  the  solution,. 
However,  as  the  number  of  elements  is  increased,  so  do  the  number  of  eigenvalues  and 
eigenfunctions  calculated. 

Now  that  the  eigenfunctions  have  been  found,  the  Franck-Condon  Factors  can  be 
calculated. 


Conclusion 

This  chapter  showed  how  the  knowledge  of  Einstein’s  Coefficients  are  vital  to 
transition  analysis  and  laser  engineering.  Einstein’s  Coefficients  can  be  derived  through 
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the  relation  of  Franck-Condon  Factors.  However,  these  factors  are  only  revealed  if  the 
molecule  is  known  intimately,  namely  through  its  eigenfunctions  which  describe  the 
nuclear  motion.  Morse  found  the  eigenfunctions  analytically  for  a  potential  energy  curve 
similar  to  a  molecule.  Unfortunately,  this  model  is  adequate  only  for  the  lowest  bound 
states  of  molecules 

The  SDO  method  showed  how  Schrodinger's  Equation  can  be  solved  numerically 
for  an  arbitrary  potential,  V(R) .  But,  if  V(R)  is  not  known  for  the  molecule,  how  do  you 
solve  for  the  molecule’s  eigenfunctions? 

The  next  chapter,  Computer  Modeling,  describes  how  the  SDO  method 
exhaustively  solves  Schrodinger's  Equation  over  and  over  again  hunting  down  that 
elusive  V(R) .  Once,  a  good  representation  of  V(R)  is  found,  the  eigenfunctions  are 
computed. 


111.  Computer  Modeling  and  Programs 


Introduction 

This  chapter  describes  how  the  SDO  method  is  implemented  on  a  computer  to 
calculate  the  Franck-Condon  Factors.  Contrary  to  the  emphasis  of  the  previous  chapter, 
the  most  important  feature  of  the  SDO  method  is  how  it  searches  and  finds  the  “best” 
potential  curve  for  the  molecule  in  question.  Thereby,  knowing  the  potential,  the 
wavefunctions  can  be  calculated  with  the  method  of  finite  elements,  which  lead  to  the 
Franck-Condon  Factors. 

Experimental  spectroscopic  data  is  central  to  the  search  for  the  “best”  potential. 
First,  this  chapter  explains  how  this  transition  data  needs  to  be  presented  to  the  computer 
model.  Next,  a  description  explains  how  this  spectroscopic  data  is  used  in  the  hunt,  that 
is,  the  search  for  the  “best”  potential.  Followed  by  a  section  which  describes  how  to  run 
the  code  which  includes  execution  strategies.  And  finally,  once  the  potential  curve  has 
been  identified,  how  to  numerically  calculate  the  Franck-Condon  Factors. 
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Background  of  Computer  Models 


The  vast  majority  of  this  research  effort  was  spent  fixing,  changing,  creating,  and 
testing  FORTRAN  code.  This  section  documents  these  efforts.  The  specific  details  of 
these  codes  is  discussed  later. 

Four  major  FORTRAN  codes  previously  written  in  1984  (2)  were  used  for  this 
project.  Unfortunately,  these  codes  survived  only  on  paper.  Therefore,  the  codes  were 
either  typed  in  by  hand  or  scanned  in  with  a  machine.  Neither  of  these  processes  worked 
all  too  well.  A  great  deal  of  time  was  spent  searching  for  bothersome  bugs  created  by 
typos  and  mistakes,  and  some  were  not  obvious.  To  add  to  the  confusion,  the  FORTRAN 
compiler  used  for  this  project  interpreted  some  of  this  old  code  differently.  Again,  more 
modifications  had  to  be  made  to  make  the  code  work  as  intended. 

Several  of  the  commercial  subroutines  used  by  these  programs  no  longer  existed 
and/or  were  missing.  These  routines  included  IMSL  Math/Library®  functions  and  a  non¬ 
linear  minimization  routine  (5).  Both  of  their  intended  functions  are  discussed  later  in 
detail. 

After  1 984,  IMSL  decided  to  completely  revamp  their  libraries.  This  meant  the 
algorithms  previously  used  did  not  exist  any  more.  In  fact,  much  of  the  code  previously 
written  had  been  specifically  adapted  to  use  these  missing  routines.  These  adaptations 
mainly  included  matrix  storage  modes  compatible  with  these  math  routines.  Therefore, 
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not  only  did  new  IMSL  routines  have  to  be  figured  out,  but  the  old  matrix  storage  modes 
within  the  code  had  to  re-written. 

Why  bother  with  these  IMSL  routines  in  the  first  place?  The  purpose  for  using 
the  IMSL  math  routines  is  simple.  This  code  requires  large  matrix  manipulations  and 
calculations.  Due  to  the  complexities  of  this  math,  such  as  diagonalizing  a  1 500x1500 
matrix,  this  commercial  product  was  selected  to  avoid  any  possible  math  problems. 

The  non-linear  minimization  routine  used  by  the  main  program  no  longer  existed 
in  digital  format  either.  This  routine,  acquired  from  Professor  Pearson  of  the  University 
of  Washington,  is  an  older  code  written  before  modem  programming  conventions  became 
standard.  This  situation  created  run-time  errors  that  were  not  initially  anticipated. 

Further,  new  code  and  subroutines  were  written  to  handle  unforeseen  complexities 
encountered  with  this  research.  The  code  previously  depended  on  a  rigid  placement  of 
the  finite  element  grid.  However,  convergence  problems  were  discovered  with  this 
approach.  Therefore,  a  new  algorithm  was  written  to  adaptively  change  the  grid  during 
code  execution.  Another  grid  problem  was  discovered  when  calculating  the  Franck- 
Condon  Factors.  These  overlap  integrals  require  the  same  grid  be  employed  for  two 
distinct  executions  of  the  code.  A  new  subroutine  was  written  to  ensure  these  two  grids 
coalesced,  even  sliding  one  of  the  grids  if  needed. 

Countless  other  small  coding  changes  were  also  made.  These  included  improved 
output  listings  and  re-dimensioning  of  variables  to  handle  larger  matrices. 
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The  rest  of  the  chapter  now  focuses  on  a  detailed  discussion  about  these  four  main 
programs.  The  remaining  sections  include  explanations  about  some  of  the  algorithms, 
“how  to  run”  the  codes,  and  execution  strategies. 
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Preparation  of  Spectroscopic  Data 


To  find  the  best  potential  surface  which  describes  the  molecular  forces, 
experimental  spectroscopic  data  must  be  supplied  to  the  model.  This  spectroscopic  data 
is  used  as  the  benchmark  for  the  model  to  decide  which  potential  “best”  matches  the 
molecule.  This  spectroscopic  data  must  be  presented  in  a  vibrational  eigenvalue  format 
for  each  electronic  state  of  interest. 

To  develop  the  vibrational  eigenvalues,  the  experimental  spectroscopic  data  can 
be  furnished  in  either  of  two  forms.  The  first  form,  and  the  most  obvious,  is  a  listing  of 
the  actual  observed  emission  or  absorption  lines  of  the  molecule.  The  other  acceptable 
way  is  to  offer  the  spectroscopic  data  in  terms  of  the  Dunham  Coefficients. 

In  1932,  Dunham  devised  a  formula  which  compactly  describes  the  energy  of  a 
rotating  vibrator  in  terms  of  coefficients  (1).  His  formula  is 


FvK  = 


(III-1) 


where  Fvk  is  the  energy  (the  eigenvalue)  of  the  rotating  vibrator,  v  is  the  vibration  level,  K 
is  the  rotational  level,  and  Yy  are  the  Dunham  coefficients. 
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Therefore,  if  the  Dunham  coefficients  are  known  for  a  given  electronic  state  of  the 
molecule,  then  the  energy  eigenvalues  can  be  calculated.  Often,  these  coefficients  are 
reported  in  journal  articles  in  lieu  of  the  actual  spectroscopic  data. 

Ostdiek  developed  a  FORTRAN  code  (2:  5-6,  44-48,  Al,  B1-B8),  called  dunham, 
which  automatically  calculates  the  energy  eigenvalues  given  these  coefficients  as  input. 
Figure  III-l  illustrates  a  sample  input  file  for  the  program  dunham. 


A  title  for  the  output  fUe^} 


>LBL=PbO  X  State 
>ROT=0 
>VIB=25 
>LVL=25 
>DEQ=0.0 
>Y00=0. 

>Y1 0=722.687 
>Y20=-3.613 


The  number  of  levels 
to  be  written  to  output  file^ 


e  maximum  number1 
of  rotation  levels 


Figure  III-l.  This  figure  illustrates  a  sample  input  file  for  the  Dunham  software  (2:47).  This 
example  only  gives  the  first  three  Dunham  coefficients  for  the  ground  electronic  state  of  PbO.  The 
accuracy  of  the  energy  eigenvalues  can  be  increased  by  supplying  more  Dunham  coefficients.  Since, 
the  rotational  fine  structure  is  not  needed,  ROT  is  always  set  to  zero.  The  output  of  dunham  are 
vibrational  eigenvalues  for  the  electronic  state  of  interest. 
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A  second  FORTRAN  code,  called  efit  (2:  6-14, 48-54,  C1-C2,  D1-D13),  takes 


observed  vibronic  transition  data  and  fits  the  data  into  an  eigenvalue  form.  This  code 
provides  an  alternative  method  over  dunham.  The  observed  spectra  from  diatomic 
molecules  often  are  the  electronic  transitions.  Within  these  transitions  is  buried  the 
vibrational  structure.  This  structure,  or  splitting  of  the  lines,  is  measured  with  a 
monochrometer  or  some  similar  device.  The  eigenvalues  can  then  be  determined  by  some 
differencing  techniques  using  these  observed  spectral  lines.  Figure  III-2  demonstrates  this 
point.  The  figure  shows  three  vibronic  transitions,  all  from  a  common  upper  vibrational 
level.  When  illustrated  such  as  this,  clearly  two  of  the  lower  vibrational  eigenvalues  can 
be  determined  by  differencing  the  observed  transitions. 
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Molecular  Transitions 


Figure  III-2.  This  figure  illustrates  vibronic  transitions  of  diatomic  molecules.  In  this  example, 
electronic  state  A  is  transitioning  to  the  ground  state  labeled  X.  The  observed  spectrum  is  indicative 
of  the  energy  difference  between  the  levels.  If  the  transitions  can  be  identified  as  to  which  are  taking 
place,  such  as  v’=3  to  v”=0,  then  the  vibrational  eigenvalues  of  each  state  can  be  determined.  For 
clarity,  only  three  transitions  are  illustrated  here.  A  simple  differencing  method  could  be  employed 
to  garner  the  first  two  eigenvalues  of  the  ground  state.  The  output  of  efit  is  a  list  of  the  vibrational 
eigenvalues  for  each  electronic  state  used  for  the  fit. 

Due  to  measurement  error  and  uncertainty  in  the  data,  several  transitions  may 
indicate  different  values  for  the  same  eigenvalue.  Obviously,  this  can’t  be  true.  The 
program  efit  uses  a  least  squares  fit  to  resolve  these  inconsistencies  between  the  data.  A 
transition  line  involving  one  energy  level  is,  therefore,  made  to  agree  with  other 
transitions  lines  involving  the  same  energy  level.  (2:  8-9)  This  correction  is  obtained  by 
minimizing  the  sum 
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s= 


( III-2) 


where  S  is  the  sum  to  be  minimized,  wy  is  the  weighting  factor,  ly  is  the  observed 
transition  line  between  state  i  and  j,  £j  and  z}  are  the  eigenvalues  of  state  i  and  j, 
respectively. 

This  program  can  accept  up  to  62,500  transition  lines,  involving  250  energy  levels 
from  10  different  electronic  states.  Figure  III-3  illustrates  for  a  sample  input  file  for  the 
program  efit. 


Associates  the  first 
level  with  label  X 


- 

>STAT01 

>STAT02: 

>LVLS01 

>SHIFTS= 

>A08X01= 

>A06X0T 

>A04X00= 

>A03X00= 

A06X02= 


Indicates  which  vibrational 
levels  should  be  included  in 
least  squares  fit 


=0  1  234567 
:-30854.15 
=22483.5 
=21675.0 
=21519.3 
=21047.8 
=20942.0 


Transition  from  state 
A  (v’=6)  to  X  (v”=2) 
sin  wavenumbers  (cm-1 


weighting  factor 
assigns  confidence 
in  observed  data 


Figure  III-3.  This  figure  illustrates  the  input  file  for  the  program  efit.  This  sample  input  file  should 
not  be  taken  literally,  a  real  input  would  include  more  transitions  than  indicated  here.  The  user 
specifies  the  energy  for  each  transition  with  it  s  assignment.  The  weighting  factor  defaults  to  1.0  if 
none  is  specified.  The  output  of  efit  is  the  energy  eigenvalues  for  each  electronic  state  used  during  the 
fit. 
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Now  that  the  eigenvalues  have  been  determined  for  each  electronic  state,  it  is  now 


time  to  find  the  “best”  potential  which  describes  this  state. 


The  Hunt 

The  program  diatom  is  responsible  for  hunting  down  the  “best”  potential  surface. 
The  “best”  potential  is  identified  by  finding  the  potential  surface  which  yield  numerical 
eigenvalues  which  are  the  closest  to  the  experimental  data  (2:  54-84,  E1-E2,  F1-F44). 

But  before  the  procedure  for  executing  diatom  is  explained,  some  background 
information  about  this  algorithm  must  first  be  addressed.  For  a  complete  description  of 
diatom ’s  input  files,  and  execution  strategies  see  the  next  section,  Executing  the  Code. 

This  background  material  includes  a  discussion  about  the  iterative  hunting 
process,  the  parameterized  potential  functions,  the  non-linear  minimization  routine,  and 
automatic  gridding. 

The  program  diatom  uses  a  non-linear  minimization  routine  to  find  the  “best” 
potential  curve.  The  minimization  routine  accomplishes  this  by  adjusting  parameters  of  a 
selected  potential  function.  Then  to  test  its  hypothesis,  Schrodinger's  Equation  is  solved 
with  the  finite  element  method  using  for  this  potential  curve  with  these  new  parameter 
values,  (see  Chapter  II).  The  difference  between  the  computed  and  experimental 
eigenvalues  then  is  used  to  provide  feedback  to  the  minimization  routine.  This  process  is 
repeated  exhaustively  until  the  optimum  potential  surface  has  been  found 
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The  following  equations  lists  the  parameterized  potential  functions  which  are 
available  to  diatom  and  the  user.  These  functions  model  the  vibrational  energy  forces 
only.  The  rotational  energy  is  not  modeled  with  this  finite  element  method  (see  Chapter 
II),  nor  do  the  following  equations  capture  the  rotational  energy.  The  parameters  of  these 
function  are  what  the  minimization  routine  adjusts  in  its  search  for  the  global  minimum. 


V(R;  k) 
V(R;  De,  a) 

V(R;  De) 

V(R;  De,  a,  P) 


1  k  (R-R,)!+  t, 

2 


D„ 


D„ 


1  _  e'°(R-R‘))2  -  i' 


( D  V2  ^ 


R^ 

vRy 


-2 


J 

6  ^ 


v  R  j 


+  T. 


+  T 


D„ 


a  -  P 


'  R  ' 


^  r; 


a 

P 


u; 


p  ^ 


+  T. 


(Harmonic) 

(Morse) 

(III-3) 

(Lennard  -  Jones) 

(Mie) 


where  V  is  the  potential  function,  R  is  the  intemuclear  separation,  k  is  the  harmonic  force 
constant,  Re  is  the  equilibrium  distance,  De  is  the  dissociation  energy,  a  is  the  Morse 
anharmonicity  term,  a  and  P  are  exponents  in  the  Mie  potential,  and  Te  measures  the 
energy  separation  between  the  minima  of  the  two  electronic  state  potential  energy  curves. 
(3:136).  The  adjustable  parameters  are  k,  De,  a,  a,  and  p. 

The  Harmonic  and  Morse  functions  are  the  only  functions  listed  in  Equation  III-3 
that  have  known  analytic  solutions.  They  were  explicitly  included  in  this  list  for  this  very 


III- 11 


reason.  Their  solutions  are  discussed  in  Chapter  II  and  are  invaluable  for  the  validation 
process.  For  more  information  about  the  validation  process  see  Chapter  IV. 

Next  on  the  list  is  the  Lennard-Jones  potential.  The  Lennard-Jones  potential  is 
cleverly  constructed  to  take  advantage  of  the  fact  that  van  der  Waals’  intermolecular 
attractive  forces  go  as  R'6.  This  attraction  is  a  “self-generating  boot-strap  force”,  called 
the  London  dispersion  force.  The  London  dispersion  force  is  caused  by  the  fluctuating 
electron  density  of  each  molecule.  In  effect,  a  temporal  dipole  moment  from  one 
molecule  induces  a  dipole  moment  on  a  neighboring  molecule,  which  in  turn  induces 
back  (4:  52-54,1 15-1 17).  The  repulsive  branch,  modeled  by  the  R'  term,  has  no 
physical  foundation.  This  repulsion  term  was  chosen  as  a  mathematical  convenience  to 
ensure  the  minimum  of  this  potential  occurs  at  Re. 

The  Mie  potential  is  formed  when  both  exponents  of  the  Lennard-Jones  potential 
are  parameterized.  Most  molecular  forces  don’t  exactly  go  as  R"6  and  R12.  By  varying 
these  exponents,  the  Mie  potential  can  be  adapted  to  match  these  actual  forces.  Both  the 
attractive  and  repulsive  branches  of  the  Mie  can  be  varied  independently,  unlike  the 
others.  Therefore,  the  Mie  potential  is  the  most  flexible  of  the  potentials  listed  in 
Equation  III-3.  The  additional  coefficients  ensure  the  minimum  of  this  potential  remains 
at  Re.  With  a  few  algebraic  substitutions,  the  Morse  and  Lennard-Jones  potentials  can  be 
shown  to  be  special  cases  of  the  Mie  potential  (2:  17). 

The  equilibrium  distance,  Re,  and  the  energy  separation,  Te,  are  not  variable 
parameters.  Within  diatom,  these  values  are  treated  as  constants.  Since  this  analysis  does 
not  include  the  rotational  fine  structure,  Re  must  be  supplied  by  the  user.  Te  must  also  be 
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supplied  by  the  user.  Te  is  defined  as  zero  for  the  ground  electronic  state.  For  details  on 
how  to  calculate  these  values  see  the  section  labeled  Executing  the  Code. 

For  the  potential  selected  by  the  user,  the  minimization  routine  (5)  adjusts  the 
potentials  parameters  in  a  systematic  manner.  First,  the  eigenvalues  are  computed  for  an 
initial  “guess”  of  the  parameters  in  question.  Then,  the  sum  of  residuals  squared  is 
calculated  by  differencing  the  experimental  eigenvalues  (the  bench  mark)  and  the  newly 
computed  eigenvalues.  The  formula  used  is 


2 

1  ^  /  \ 

Slim  —  ^  1  experimental  ^calculated  )  (III-4) 

Z  j  =  l 

where  wj  is  the  weight  factor,  Experimental  is  the  observed  eigenvalue,  E'caiCuiated  is  the 
computed  solution,  and  n  is  the  maximum  number  eigenvalues  furnished  by  the  user.  The 
minimization  routine  then  adjusts  the  parameters  of  the  chosen  potential  (hopefully  in  the 
direction  to  the  global  optimum),  the  eigenvalues  are  then  re-computed,  the  sum  of 
residuals  squared  is  recomputed,  and  so  on.  This  process  is  exhaustively  repeated  until 
the  sum  of  the  residuals  squared  has  been  minimized,  thereby,  the  “best”  solution  has 
been  isolated.  See  Figure  III-4. 
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Figure  III-4.  This  figure  illustrates  schematically  the  iterative  hunting  process.  From  an  initial 
starting  place  given  by  the  user,  the  model  computes  the  vibrational  eigenvalues  using  the  method  of 
finite  elements.  The  sum  of  the  residuals  squared  is  then  computed  (See  Equation  III-4)  and  the  non¬ 
linear  minimization  routine  then  determines  which  direction  must  be  made  to  find  the  minimum.  The 
minimization  routine  then  adjusts  the  parameters  of  the  potential  accordingly,  and  reseeds  the 
process.  The  process  is  complete  when  the  global  minimum  has  been  found,  that  is,  the  best  potential 
which  yields  the  eigenvalues  closest  to  experimental  results. 


This  non-linear  minimization  routine  searches  for  the  global  minimum  of  an  N- 
dimensional  parameter  space  by  making  a  series  of  gradient  steps,  random  direction  steps, 
average  directional  steps,  and  random  jumps  within  the  parameter  space.  This  program 
uses  the  method  of  steepest  descents  for  its  gradient  maneuvers.  If  the  routine  hits  a 
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boundary  of  the  parameter  space  it  will  “slide”  along  the  boundary  until  improvement 
ceases.  The  random  jumps  ensure  the  routine  does  not  fall  into  a  local  minimum,  thereby, 
missing  the  global  optimum. 

The  parameter  space  in  which  this  minimization  routine  searches  for  the  global 
optimum  is  graphically  shown  in  Figure  III-5.  This  surface  was  generated  by  calculating 
the  sum  of  the  residuals  squared  for  over  1000  points.  The  analytical  eigenvalues  of  the 
Morse  Oscillator  were  provided  as  the  benchmark  for  a  dissociation  energy  (De)  of 
15000.0,  and  an  anharmonicity  term  (a)  of  6.0.  Then  the  sum  of  the  residuals  squared 
was  calculated  as  a  was  varied  from  1 .0  to  12.0,  and  De  was  varied  from  6000.0  to 
24000.0  .  The  height  of  the  surface  in  Figure  III-5  is  the  log  of  the  sum  of  residuals 
squared  at  each  (a,  De)  coordinate  pair.  This  height  was  then  truncated  so  the  lowest 
value  was  6.  This  truncation  was  done  in  order  to  see  more  resolution  in  the  plot. 

The  graph  shows  some  very  important  features  of  this  parameter  space.  The  most 
important  feature  is  that  there  is  only  one  minimum.  And  this  minimum  is  extremely 
pronounced.  Another  feature  is  that  the  surface  which  the  minimization  routine  moves 
upon  appears  smooth,  for  the  most  part.  A  gradient  search  method  should  be  sufficient  to 
find  the  global  minimum  for  a  surface  such  as  this.  The  two  triangular  bumps  are  real 
and  were  calculated  independent  of  each  other.  The  origin  of  these  bumps  is  not  currently 
understood. 
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Figure  III-5.  This  figure  illustrates  the  parameter  space  for  where  the  hunt  occurs.  The 
minimization  routines  roams  around  on  this  surface  until  it  finds  the  global  optimum.  The  analytical 
solutions  of  the  Morse  potential  where  a= 6.0  and  D,.=15000.0  were  provided  as  the  benchmark.  This 
benchmark  included  all  28  bound  states.  The  height  of  this  surface  represents  the  order  of 
magnitude  (log  scale)  of  the  sum  of  residuals  squared  (i.e.  a  measure  of  how  close  the  computed 
eigenvalues  are  to  the  benchmark).  The  sum  of  the  residuals  squared  was  computed  for  a  from  1.0  to 
12.0  (the  x-axis),  and  for  De  from  6000.0  to  24000.0  (the  y-axis)  for  300  grid  elements.  The  bottom  of 
this  well’s  order  of  magnitude  is  actually  101,  it  was  truncated  at  106  for  illustrative  purposes. 
Another  important  point,  because  this  is  an  order  of  magnitude  graph,  subtle  changes  in  value  may 
be  washed  out.  The  two  triangular  bumps  were  calculated  independent  of  each  other. 


The  math  routines  which  diagonalizes  H  matrix  (see  Chapter  II)  are  two  IMSL© 
Math/Library™  FORTRAN  subroutines  (6:  307-311).  The  first  routine,  DEVLSF, 
computes  the  eigenvalues  of  a  real  symmetric  matrix  by  an  orthogonal  similarity 
transformation  to  an  equivalent  symmetric  tridiagonal  matrix,  then  it  performs  an  implicit 
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QL  algorithm  to  compute  the  eigenvalues  of  this  matrix.  The  second  routine,  DEVCSF, 
does  the  same  thing,  but  also  returns  the  eigenfunctions.  These  eigenfunctions  are 
normalized  such  that  the  infinity  norm  of  each  vector  is  one. 

The  important  point  here  is,  while  the  minimization  routine  is  hunting  for  the  best 
potential,  DEVLSF  is  the  routine  that  is  used  to  compute  the  eigenvalues.  This  saves  time 
because  DEVLSF  works  faster  than  DEVCSF.  Once  the  best  potential  has  been  found, 
then  DEVCSF  is  called  to  calculate  the  eigenfunctions.  Therefore,  DEVCSF  is  only 
called  once  at  the  very  end  of  execution. 

In  the  original  version  of  diatom,  the  user  had  complete  control  over  the  finite 
element  grid.  The  user  specified  where  the  two  grid  endpoints  were  located,  and  the 
number  of  elements  to  be  used.  Once  the  grid  was  defined,  however,  this  grid  was  used 
for  all  calculations,  no  matter  how  the  minimization  might  adjust  the  parameters  of  the 
potential.  Refer  to  Figure  II-7. 

The  convergence  of  the  finite  element  method,  as  it  turned  out,  is  very  sensitive  to 
how  the  grid  was  laid  down.  See  Chapter  IV  for  a  discussion  on  convergence.  For 
example,  if  the  endpoints  of  the  grid  were  positioned  so  that  only  a  small  part  of  the 
potential  was  included  for  the  calculation,  the  solution  would  not  to  converge.  This 
situation  is  easily  obtained.  Recall  the  global  minimization  routine  adjusts  the  parameters 
of  the  potential  function.  And  these  adjustments  may  be  systematic  or  even  random. 
Therefore,  the  potential  curve  may  swing  completely  off  a  previously  defined  grid. 

Therefore,  the  program  diatom  has  been  modified  in  this  work  to  automatically 
adjust  the  grid  during  execution  for  the  user.  The  endpoints  of  the  grid  are  adjusted  so 
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that  the  potential  curve,  and  the  bound  state  wavefunctions  do  not  get  truncated  by  the 
grid  endpoints  as  the  minimization  adjusts  the  parameters  of  the  potential  function.  This 
smart  grid  algorithm  automatically  “walks  the  endpoints  in  and  out”  every  time  the 
potential  function’s  parameters  are  changed  by  the  minimization  routine. 

The  endpoints  are  adjusted  so  that  the  value  of  the  potential  satisfies  certain 
conditions.  The  value  of  the  repulsive  branch  at  the  left  endpoint  must  fall  between  10 
times  and  1000  times  De.  The  value  of  the  attractive  branch  at  the  right  gridpoint  must 
fall  between  -10"1  and  -10'6  times  De.  These  values  were  chosen  from  experience  and 
experiment.  Other  than  the  fact  these  limits  work,  they  were  chosen  arbitrarily. 


Executing  the  Code 

This  section  does  two  things.  First,  the  input  files  are  explained  in  detail.  Second, 
execution  strategies  are  discussed. 

The  program  diatom  requires  two  input  files  for  execution.  One  input  file  defines 
the  bench  mark  for  the  best  potential,  that  is  the  experimental  vibrational  eigenvalues  for 
the  electronic  state  of  interest.  The  second  input  file  is  a  master  control  file  which 
manages  the  execution. 
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Figure  HI-6  shows  a  sample  eigenvalue  energy  input  file. 


diatom_energy  .input 


>LBL=MORSE:  De=15000  cm-1,  a  =  6.0 
>V00=-14484. 8847577293 
1.654273 1880 
7886467 
'33041054 
7.9628195640 


Title  to  be  Included 
in  Output 


VU7 

>V05 


TT71 

-9828.7323350227 


>V23=-5 18.58361 3. 
>V24=-343.35312 
>V25=-204. 122644 
>V26=- 100.892 1596549 
>V27=-33. 66 1675 1135 
>V28=-2.43 11 905722 


Vibrational  Energy 
Eigenvalue  at  v=5 


Figure  III-6.  This  figure  illustrates  the  energy  input  file  for  diatom  (“diatom_energy.input”).  This 
file  defines  the  benchmark  for  the  best  potential.  The  model  will  attempt  to  find  a  potential  whose 
eigenvalues  match  most  closely  to  these.  As  an  example,  this  list  includes  the  28  analytical 
eigenvalues  for  a  Morse  Oscillator  where  a=6.0  and  De=15000.0.  The  format  for  each  energy  input  is 
“>Vxx”,  where  xx  is  the  state  label.  Labels  less  than  10  still  must  have  digits,  such  as  “>V05”. 
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Figure  HI-7  shows  the  master  control  input  file. 


(jTof  writes)—--^^ 

/#ofgradient 
V^jries 

\Random  Jumpy 
(Run  WaviT 


(^Reduced  Mass^> 
(Qonstant  #lRe^. 


Constant  #2  Te. 


diatom_control.input 

>GP=1- 
>BG=-4.0 
>EN=8.0 
>NE=300 
>SG=0 
>GA=1 
>RM=0 
>NT=40 
>NW=0 
>NR=1 
>NG=1 
>NA=1 
>NJ=1 
>FR=1 
>FP=1 
>RW=1 
>HB=1. 

>MU=1.0 
>L  1=3.0 
>P  1=4.0. 

>U  1=5.0 
>C  1=2.0 
>C2=0.0 


Figure  III-7.  This  figure  illustrates  the  master  control  input  file  for  diatom  (“diatom  control.input”). 
This  file’s  syntax  is  “>AA=“  where  AA  is  a  two  letter  code.  AH  of  the  possible  codes  are  exemplified 
in  this  sample.  This  input  files  consists  of  three  main  groups.  The  first  group  controls  the  grid,  the 
second  group  controls  the  minimization  routine,  and  the  third  group  controls  the  potential  and  its 
associated  parameters  and  constants.  The  two  letter  code,  RM,  turns  the  non-linear  minimization 
routine  on  and  off. 
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The  two  letter  code  GP  selects  the  generalized  potential  as  shown  Table  ID-1 . 


Table  III-l.  This  table  shows  how  to  select  the  GP  flag. 


Potential  Function 

GP  value 

Simple  Harmonic  Oscillator 

1 

Morse  Anharmonic  Oscillator 

2 

Lennard- Jones  Potential 

3 

Mie  Potential 

4 

Each  potential  function  has  a  number  of  parameters  and  constants  (refer  to 
Equation  III-3).  The  user  has  the  ability  to  define  a  search  space  for  each  parameter.  That 
is,  define  the  upper  and  lower  limits  for  each  parameter,  and  define  its  initial  starting 
value.  These  values  define  the  parameter  space  where  the  minimization  routine  hunts  for 
the  best  potential. 


Table  III-2.  This  table  shows  the  relationship  between  the  parameter  and  constant  variables,  and  the 
numbers  they  are  identified  by  in  the  control  input  file. 
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For  example,  the  lower  limit,  upper  limit,  and  initial  value  of  a  (Morse’s 
anharmonicity  term)  would  be  denoted  by  the  codes  L2,  P2,  and  U2,  respectively.  See 
Figure  HI-7  for  the  example  input  file. 

Optimally,  the  user  wants  to  initialize  the  parameters  as  close  to  the  expected 
values  as  possible.  This  preemptive  measure  will  help  the  search  routine  find  the  “best” 
potential  sooner.  Also,  set  the  lower  and  upper  limits  as  tight  as  possible.  This  will  limit 
the  search  volume,  and  again  speed  things  up.  Later,  Chapter  IV  will  show  that  solutions 
of  acceptable  accuracy  require  a  lot  of  elements,  which  translates  into  lots  of  time.  Some 
quick  estimates  of  the  solution  can  be  extremely  beneficial. 

Since  this  analysis  does  not  include  the  rotational  fine  structure,  Re  must  be 
supplied  by  the  user.  The  equilibrium  distance.  Re,  can  be  computed  by  using  the 
rotational  spectral  constant  Be.  Te  may  be  approximated  by  finding  the  difference 
between  the  ground  vibrational  eigenvalue  of  the  upper  electronic  state  and  the  ground 
vibrational  eigenvalue  of  the  ground  state  However,  a  preferable  solution  is  to  calculate 

T' by  subtrac,in8  ground  vibrational  dgenvalue  of each 

electronic  state,  and  then  subtracting  these  two  values.  This  method  assumes  the 
potential  energy  curve  behaves  like  a  Morse  Oscillator  at  the  bottom  of  the  well,  which  in 
most  cases  is  a  good  approximation.  Again,  coe  and  ©eXe  are  vibrational  constants. 

The  grid  is  controlled  by  five  settings.  The  settings,  BG  and  EN,  position  the  left 
and  right  grid  endpoints  respectively.  These  values  have  dimensions  of  length.  The 
setting  NE  carves  up  the  grid  into  the  number  of  elements  specified,  all  of  equal  size. 
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The  setting  SG,  the  smart  grid  algorithm,  has  three  settings.  Setting  this  parameter  to  0 
turns  it  off,  1  and  2  turns  it  on.  With  SG  turned  on,  the  smart  grid  walks  the  grid 
endpoints  in  and  out  as  needed  to  keep  the  grid  space  defined  properly  for  the  potential 
curve  in  question.  With  SG  set  to  1,  the  initial  stepsize  is  kept  constant,  therefore,  the 
number  of  elements  vary  as  the  endpoints  are  adjusted.  With  SG  set  to  2,  the  number  of 
elements  is  kept  constant,  therefore,  the  stepsize  varies  as  the  endpoints  are  adjusted.  The 
setting  GA,  the  grid  alignment  algorithm,  co-aligns  grids  of  separate  executions.  Setting 
it  to  0  turns  it  off,  and  setting  it  to  1  turns  it  on.  This  algorithm  records  previously 
defined  grids  by  writing  data  to  a  special  output  file.  This  information  is  then  used  to 
align  the  grids.  The  grid  alignment  algorithm  must  be  turned  on  if  the  eigenfunctions  are 
going  to  be  used  for  Franck-Condon  calculations.  Note,  turning  the  grid  alignment  on 
overrides  the  smart  grid  assignment  of  2,  switching  it  to  1 . 

The  non-linear  minimization  routine  has  7  settings  to  control  its  execution.  The 
flag  RM  turns  the  minimization  routine  on  and  off  with  1  and  0,  respectively.  If  the 
minimization  routine  is  off,  the  code  will  calculate  the  eigenvalues  for  the  initial 
parameter  settings  only.  When  on,  the  minimization  routine  will  make  NT  number  of 
attempts  to  find  the  global  minimum.  The  flags  NR,  NG,  NA,  NJ  control  the  relative 
number  out  of  NT  tries  each  of  these  strategies  will  be  attempted.  The  setting  NR 
controls  the  number  of  times  it  tries  a  random  direction  looking  for  the  minimum,  NG  is 
the  number  of  times  gradient  steps  should  be  taken,  NA  is  the  number  of  average 
directions  it  should  take,  and  NJ  is  the  number  of  times  it  should  take  a  random  jump  to 
some  other  point  in  the  parameter  space.  NW  controls  the  number  of  times  it  reports 
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diagnostic  information  to  the  screen.  It  will  never  provide  diagnostics  if  it  is  set  to  0,  but 
will  report  back  every  other  step  if  the  setting  is  2,  etc. 

The  designators  HB  and  MU  set  the  values  of  Planck’s  Constant  over  2n  (h  bar) 
and  the  reduced  mass  of  the  molecule.  Units  are  important.  The  grid  endpoints,  BG  and 
EN,  HB,  MU,  and  the  energy  eigenvalues  all  need  to  use  the  same  unit  system.  Atomic 
units  are  recommended,  because  atomic  units  allow  HB  to  be  set  to  one,  thereby, 
reducing  numerical  error  introduced  by  using  very  small  numbers.  As  an  example,  in  SI 
units  Planck’s  constant  is  on  the  order  of  10'34.  Using  such  a  small  number  could  lead  to 
significant  numerical  error,  especially  when  using  any  algorithm  to  diagonalize  matrices 
such  as  DEVLSF. 

The  last  settings,  RW,  FP,  and  FR  are  option  flags  which  turn  on  and  off  with 
settings  of  1  and  0,  respectively.  RW  instructs  the  diagonalization  math  routine  to  return 
both  the  eigenfunctions  and  eigenvalues.  This  setting  saves  time  for  the  user  when  turned 
off  if  the  user  is  not  interested  in  the  eigenfunctions.  FP  and  FR  instruct  the  model  to 
save  a  potential  and  residual  output  file.  The  potential  file  contains  1000  coordinate  pairs 
of  position  versus  the  value  of  the  “best  found”  potential  at  each  position.  The  residual 
file  contains  the  listing  of  residuals  for  each  experimental  eigenvalue  compared  to  its 
associated  “best  found”  numerical  eigenvalue. 

A  good  strategy  for  running  diatom  is  to  start  simple.  Rims  with  lots  of  elements 
are  expensive.  The  global  minimization  routine  may  solve  Schrodinger's  Equation 
hundreds  of  times  depending  on  the  execution  and  its  input.  Therefore,  starting  with  a 
large  search  volume  and  a  lot  of  elements  could  be  costly.  Start  with  not  as  many 
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elements,  like  a  150  or  200.  As  the  minimization  routine  starts  to  hunt  down  the  best 
potential,  winnow  the  search  volume.  Also,  run  with  the  smart  grid  algorithm  on.  An 
initial  setting  of  2  will  help  control  the  number  of  elements  used,  and,  therefore  control 
time.  The  hunt  for  the  best  solution  may  take  5  to  10  runs,  each  run  a  refinement  over  the 
previous. 

From  experience,  several  of  the  minimization  options  should  be  turned  off.  These 
options  include  NR  (random  directions)  and  NA  (average  direction).  Both  of  these 
options  redirect  the  minimization  routine  away  from  the  direction  to  the  global  optimum. 
Refer  to  Figure  III-5  for  justification  of  this  recommendation.  Figure  III-5  suggests  the 
parameter  space  is  smooth  and  only  one  minimum  exists.  With  this  figure  in  mind,  only 
gradient  steps  should  be  required  to  find  this  minimum.  Random  jumps  to  a  new  location 
in  the  parameter  space  may  help  avoid  identifying  the  local  minima  as  the  global  minima 
if  they  exist.  The  minimization  routine  identifies  a  minimum  by  comparing  previous 
values.  It  could  very  well  find  a  local  minimum  created  by  numerical  noise.  Therefore, 
the  recommendations  are  to  set  NR  and  NA  to  zero,  NG  to  5,  and  NJ  to  1 . 

The  program  diatom  produces  several  output  files.  The  first  one  is  diatom,  output 
which  is  the  main  output  file.  Another  is  input  Jor  Jcf. output,  this  file  contains  the 
wavefunctions  and  S  matrix  required  for  the  Franck-Condon  Factor  calculations.  The  last 
two  output  files  are  plot potential,  output  and  plot  residual,  output.  These  files  were 
described  earlier. 
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Once  the  best  potential  has  been  found  for  both  of  the  electronic  states  of  interest, 


the  Franck-Condon  Factors  can  be  calculated. 


Franck-Condon  Factor  Calculations 

This  section  describes  how  to  calculate  the  Franck-Condon  Factors  with  the 
program  called  fcf  (2:  36-40,  66-70,  G1-G2,  H1-H10).  This  program  uses  the  numerical 
wavefimctions  obtained  by  executing  diatom  and  its  finite  element  method  algorithm. 

To  calculate  the  Franck-Condon  Factors,  the  overlap  integral  for  every 
combination  of  upper  electronic  and  lower  electronic  vibrational  states  must  be 
computed.  The  eigenvectors  garnered  through  the  finite  element  method  contain 
information  about  the  wavefunction’s  value  and  slope  at  each  gridpoint.  However,  this 
data  just  runs  in  successive  order,  and  is  not  explicitly  labeled  for  each  gridpoint. 
Therefore,  information  about  the  grid  is  required  to  compute  the  overlap  integral.  This 
information  is  found  in  the  S  matrix  as  developed  in  Chapter  II. 

The  Franck-Condon  Factors  are  then  calculated  using  the  following  formula: 


where  qvV,  is  the  Franck-Condon  Factor,  v'  is  one  of  the  wavefunctions  for  the  upper 
electronic  state,  v"  is  one  of  the  wavefunctions  for  the  lower  electronic  state,  and  S  is  the 
matrix  which  describes  the  grid  space. 

After  a  successful  completion  of  diatom  (with  RW=1)  the  program  automatically 
saves  a  file  which  contains  the  first  25  wavefunctions  and  the  S  matrix  in  the  file 
input  JorJcf  output.  Two  of  these  output  files  are  required  for  the  Franck-Condon 
calculation,  one  for  each  electronic  state.  For  this  algorithm  to  work,  both  files  must  have 
the  same  S  matrix.  The  automatic  grid  alignment  algorithm  (GA=1)  must  be  invoked  for 
both  diatom  runs  to  ensure  they  both  have  the  same  grid  spacing.  The  S  matrix,  which 
depends  on  the  number  of  elements  and  grid  spacing,  is  adjusted  to  work  for  the  largest 
wavefunction. 

This  grid  alignment  algorithm,  new  to  diatom,  does  two  things.  One,  it  makes 
sure  both  runs  have  the  same  stepsize.  And  two,  if  the  gridpoints  are  not  co-aligned,  this 
algorithm  slides  one  of  meshes  over  to  coalesce  with  the  other.  In  that  way,  both  output 
files  are  ensured  to  have  the  same  S  matrix. 

The  FORTRAN  program ,fcf  takes  these  two  output  files  (renamed  by  the  user  to 
fcf  low.  input  and  fcf  high,  input)  and  calculates  the  Franck-Condon  factors.  This 
program  has  been  modified  to  handle  wavefunctions  derived  from  grids  with  different 
endpoints.  The  algorithm  builds  a  25x25  table  of  factors  by  looping  over  Equation  III-5. 
The  vector-matrix  multiplication  is  carried  out  with  two  IMSL®  Math/Library™ 
FORTRAN  subroutines  (6:  993-995).  The  first  one  is  DMURRV which  multiplies  a  real 
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rectangular  matrix  by  a  vector.  The  other  is  DMURBV  which  multiplies  a  real  band 
matrix  in  band  storage  mode  by  a  real  vector. 

The  final  output  is  the  table  of  Franck-Condon  Factors  called  fcf  output. 


Conclusion 

This  chapter  explained  how  starting  with  little  information,  namely  some 
spectroscopic  data,  the  Franck-Condon  Factors  can  be  calculated.  The  most  important 
feature  of  this  method,  however,  was  not  the  calculation  of  the  molecule’s  wavefunction 
which  give  you  these  factors.  But.  how  it  searched  for  and  found  the  “best”  potential 
curve  for  the  molecule  in  question. 
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IV.  Validation 


Introduction 

This  chapter  demonstrates  that  Schrodinger's  Equation  can  be  solved  using  the 
numerical  method  of  finite  elements.  This  validation  includes  comparisons  between  the 
analytical  and  numerical  solutions  of  the  Harmonic  and  Morse  potentials.  These  sections 
are  then  followed  by  a  discussion  about  the  convergence  of  this  method  and  its  associated 
cost.  And  finally,  the  eigenvalues  from  a  numerical  H2  molecule  potential  are  used  to 
demonstrate  how  the  model  can  search  and  find  an  unknown  potential. 

Validation  is  an  important  effort  for  any  model.  No  model  should  be  trusted  until 
it  has  the  proven  ability  to  predict  known  solutions  accurately.  Since  the  method  of  finite 
elements  has  not  been  used  a  lot  in  quantum  mechanics,  special  emphasis  is  taken  in  this 
chapter  to  demonstrate  the  validity  of  this  technique. 


Simple  Harmonic  Oscillator 

The  harmonic  oscillator  is  one  of  a  few  potentials  known  to  have  analytical 
solutions.  The  use  of  harmonic  oscillator,  therefore,  is  an  obvious  choice  for  validating 
any  numerical  solution  of  Schrodinger's  Equation.  This  section  details  comparisons 
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between  the  analytical  eigenvalues  and  functions,  and  their  numerical  equivalents 
computed  with  the  method  of  finite  elements. 

The  harmonic  oscillator,  as  modeled  in  Equation  III-3,  was  numerically  evaluated 
using  diatom.  The  force  constant  (k)  was  chosen  to  be  4.0,  and  the  reduced  mass  (p), 
h!2n,  and  the  radius  of  equilibrium  (Re)  were  all  set  to  1 .0  .  The  minimization  routine 
was  not  used. 

A  finite  element  grid  was  laid  down  over  this  potential  energy  curve  with  the  left 
endpoint  at  -5.0  and  the  right  endpoint  at  7.0  .  This  grid,  therefore,  is  symmetric  about 
Re.  A  solution  for  300  elements  was  arbitrarily  chosen  for  this  demonstration. 

With  the  parameter  choices  listed  earlier,  the  analytical  eigenvalues  conveniently 
become  1 .0,  3.0,  5.0,  and  so  on  using  Equation  11-14.  These  values  were  then  input  to  the 
diatom  as  the  bench  mark  for  residual  calculations. 

The  residuals  for  the  first  25  states  yielded  values  on  the  order  of  1  O'4.  That  is,  the 
difference  between  the  analytical  and  numerical  eigenvalues.  The  agreement 
demonstrated  by  the  relative  errors  proved  just  as  good.  They  ranged  in  value  from 
4.0xl0‘2  to  8.0x1  O'4  percent  as  shown  in  Table  IV- 1.  This  excellent  agreement  was  found 
to  be  ubiquitous  for  all  values  of  force  constants  and  reduced  masses  tested. 
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Table  IV-1.  This  table  demonstrates  the  accuracy  of  the  finite  element  method  when  solving 
Schrodinger's  Equation  for  the  Simple  Harmonic  Oscillator.  For  this  test,  k=4.0,  h/2rt=l  .0,  and 
p=1.0.  A  residual  is  defined  as  the  difference  between  the  expected  and  calculated  eigenvalues. 


Eigenstate  Level 

Calculated  Eigenvalue 

Residual 

Relative  Percent  Error 
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T 

1 

1 

1 

0.751E-03 

Sum  of  Residuals  Squared 

1 

0.2245 169E-05  1 

The  sum  of  the  residuals  squared  was  on  the  order  of  10'6.  For  this  comparison, 
the  minimization  routine  had  been  turned  off.  However,  if  a  search  had  been  conducted, 
this  value  would  have  identified  the  global  minimum. 

The  eigenfunctions  also  exhibited  excellent  agreement  with  the  analytical 
solutions.  These  eigenfunctions  were  obtained  using  the  following  method.  The 
input  Jor 'Jef  output  file  contains  a  the  listing  for  each  eigenfunction.  Eigenvectors  for 
states  0, 1,6,  and  9  were  then  selected  from  this  file.  Two  even  states  and  two  odd  states 
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were  chosen  for  their  symmetry  and  asymmetry  properties,  respectively.  Each  vector  is 
represented  by  a  long  column  of  numbers.  This  column  is  a  sequential  listing  of  pairs  of 
numbers,  each  pair  is  the  wavefunction’s  value  and  slope  at  each  gridpoint,  in  successive 
order.  For  this  analysis,  only  the  wavefunction’s  value  at  each  gridpoint  was  kept.  The 
slopes  were  discarded. 

These  four  eigenfunctions  were  then  tested  for  proper  normalization  and 

orthogonality.  All  integrations  were  conducted  using  Simpson’s  rule.  Normalization  was 

tested  by  integrating  each  wavefunction  squared  over  the  finite  element  grid  space.  All 

four  wavefunctions  yielded  a  normalization  factor  of  1 .0,  even  though  only  7  significant 

digits  were  used  to  represent  each  point.  Because  these  wavefunctions  are  from  the  same 

electronic  state,  they  are  expected  to  be  also  orthogonal  to  each  other.  The  product  of 

each  vector  with  another  vector,  including  all  permutations,  was  then  integrated.  This 

© 

integration  resulted  in  no  value  greater  than  2.92x10  .  When  considering  precision,  these 
numbers  all  can  be  interpreted  as  zero,  thereby,  demonstrating  orthogonality.  Therefore, 
all  four  eigenfunctions  tested  were  properly  normalized  and  orthogonal. 

These  four  orthonormal  eigenvectors  were  then  plotted  with  Mathematica®.  For 
comparison  purposes,  the  analytical  wavefunctions  were  also  plotted  on  the  same  graph. 
Because  these  two  wavefunctions  coalesced  so  closely  together,  the  analytical  solutions 
were  multiplied  by  -1  to  separate  the  curves.  Otherwise,  the  curves,  to  within  plotting 
quality,  would  lay  right  on  top  of  each  other.  Their  excellent  agreement  can  be  observed 
in 
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Figure  IV-l  and  in  Figure  IV-2.  Notice  how  the  amplitudes  are  the  same  for  the  numerical 
and  analytical  functions.  Also  notice  their  shapes  and  where  they  cross  the  horizontal 
axis  are  identical. 


Figure  IV-l.  This  figure  shows  the  excellent  agreement  between  the  numerical  eigenfunctions  and 
analytical  eigenfunctions  of  the  Simple  Harmonic  Oscillator.  Each  solid  curve  is  the  numerical 
solution,  and  each  dashed  curve  is  the  analytical  solution.  For  illustrative  purposes,  the  analytical 
solution  has  been  reflected  about  the  horizontal  R-axis.  The  top  curves  are  for  the  ground  state 
(n=0),  and  the  bottom  curves  are  for  the  first  excited  state  (n=l).  Note  the  symmetry  about  R^l.O  . 
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Figure  IV-2.  This  Figure  also  shows  the  excellent  agreement  between  the  numerical  eigenfunctions 
and  analytical  eigenfunctions  for  the  Simple  Harmonic  Oscillator.  Each  solid  curve  is  the  numerical 
solution,  and  each  dashed  curve  is  the  analytical  solution.  For  illustrative  purposes,  the  analytical 
solution  has  been  reflected  about  the  horizontal  R-axis.  The  top  curves  are  for  the  state  n=6,  and  the 
bottom  curves  are  for  the  state  n=9.  Note  the  symmetry  about  Re=1.0  . 
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The  Morse  Anharmonic  Oscillator 


The  Morse  potential  is  another  function  with  known  analytical  solutions.  And  as 
an  added  benefit,  because  the  Morse  function  is  such  an  excellent  approximation  of  the 
molecular  forces,  this  function  provides  a  more  rigorous  test  than  the  Harmonic 
Oscillator. 

The  validation  process  for  the  Morse  potential  proved  to  be  more  challenging  than 
that  encountered  with  the  Harmonic  Oscillator.  In  fact,  during  the  validation  and  analysis 
of  this  function,  the  requirement  for  a  subroutine  to  automatically  correct  grid  placement 
was  realized.  The  reason  for  the  need  for  an  automatic  grid  involves  convergence  issues 
which  is  discussed  in  the  next  section. 

The  Morse  Oscillator’s  analytical  eigenvalues  and  functions  were  calculated  using 
Equation  11-16.  A  dissociation  energy  (De)  of  15000.0  and  an  anharmonicity  term  (a)  of 
6.0  was  chosen  for  this  demonstration.  Again,  h/27t  and  the  reduced  mass,  p,  were  set  to 
1 .0  .  These  parameters  bear  no  resemblance  to  any  real  molecule.  But,  the  choice  of 
these  two  parameters  yields  28  bound  states,  which  is  reasonable  for  a  molecule. 

A  finite  element  grid  of  300  elements  was  arbitrarily  chosen.  The  grid  endpoints 
were  located  at  1 .5  and  2.5,  and  the  radius  of  equilibrium  (Re)  was  chosen  to  be  located  at 
2.0. 

The  residuals,  the  difference  between  the  numerical  and  analytical  eigenvalues, 
ranged  in  value  from  10'6  to  10'3,  except  for  the  last  two  bound  states  which  faired  a  little 
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worse.  Looking  at  the  relative  errors  in  Table  IV-2,  these  values  ranged  fromlO'8  to  1CT4 
percent,  again  except  for  the  last  two  states.  The  states  very  near  the  dissociation  limit 
always  proved  more  difficult  to  calculate.  These  states  are  on  the  cusp  of  being  free 
states  which  means  the  grid  needs  to  have  infinite  extent  to  model  these  wavefunctions. 


Table  IV-2.  This  table  verifies  the  accuracy  of  the  numerical  solution  of  Schrddinger's  Equation  for 
the  Morse  Potential.  For  brevity,  only  selected  levels  of  the  28  bound  states  are  shown. 


Eigenstate  Analytical 

Level  Eigenvalue 
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Calculated  Residual  Relative 
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The  numerical  wavefunctions  were  then  developed  for  states  0,  1,6,  and  9  in  a 
similar  fashion  as  described  for  the  harmonic  wavefunctions.  These  four  Morse 
wavefunctions  were  then  tested  for  normalization  and  orthogonality. 
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The  integration,  again  using  Simpson’s  Rule,  yielded  normalization  factors  of 
0.990000,  0.990001,  0.990015,  and  0.989975  for  states  0, 1,6,  and  9,  respectively.  These 
results  can  be  attributed  to  the  fact  that  only  around  130  out  of  the  300  gridpoints  are 
significant  to  the  integral.  Even  though  a  grid  was  required  with  endpoints  at  1.5  to  2.5 
for  convergence,  these  wavefunctions  only  exist  between  the  points  1 .85  and  2.29  .  The 
test  of  orthogonality  yielded  results  not  greater  than  2.34x1 0‘7,  except  for  the  integrals 
involving  state  6.  These  values  averaged  around  2.33x1 0'3.  Again,  these  numbers, 
except  for  state  6,  can  be  considered  zero  due  to  precision.  The  integrals  involving  state  6 
can  be  argued  to  almost  vanish.  Therefore,  considering  the  approximations  made  with 
Simpson’s  rules  and  the  sparse  data,  these  functions  are  orthonormal. 

Figure  IV-3  and  Figure  IV-4  show  the  excellent  agreement  with  the  analytical 
solutions.  Again,  the  analytical  solutions  were  multiplied  by  -1  to  separate  the  curves. 
Otherwise,  they  would  literally  lay  right  on  top  of  each  other.  Notice  the  amplitudes  of 
the  analytic  and  numerical  solutions  are  virtually  identical.  The  wavefunctions  also  cross 
the  horizontal  exactly  at  the  same  point. 
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Figure  IV-3. .  This  figure  shows  the  excellent  agreement  between  the  numerical  eigenfunctions  and 
analytical  eigenfunctions  of  the  Morse  Anharmonic  Oscillator.  Each  solid  curve  is  the  numerical 
solution,  and  each  dashed  curve  is  the  analytical  solution.  For  illustrative  purposes,  the  analytical 
solution  has  been  reflected  about  the  horizontal  R-axis.  The  top  curves  are  for  the  ground  state 
(n=0),  and  the  bottom  curves  are  for  the  first  excited  state  (n=l).  Note  the  asymmetry  about  Rt=2.0 . 
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Figure  IV-4.  This  figure  also  shows  the  excellent  agreement  between  the  numerical  eigenfunctions 
and  analytical  eigenfunctions  for  the  Morse  Anharmonic  Oscillator.  Each  solid  curve  is  the 
numerical  solution,  and  each  dashed  curve  is  the  analytical  solution.  For  illustrative  purposes,  the 
analytical  solution  has  been  reflected  about  the  horizontal  R-axis.  The  top  curves  are  for  the  state 
n=6,  and  the  bottom  curves  are  for  the  state  n=9.  Note  the  asymmetry  about  Re=2.0 . 
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Convergence 


This  section  details  convergence  testing  done  for  the  Harmonic  and  Morse 
potentials.  Convergence,  in  the  case  of  the  method  of  finite  elements,  means  that  the 
solution  should  approach  the  “correct”  answer  as  more  elements  are  added  to  the  grid, 
thus  reducing  the  stepsize.  However,  there  is  a  penalty  for  adding  elements  to  the  grid. 
This  penalty  is  the  cost  of  time.  In  fact,  as  shown  later,  the  time  required  for  the 
calculations  increases  cubicly  as  the  number  of  elements  are  increased. 

The  solutions  for  the  Simple  Harmonic  Oscillator  were  first  tested  for 
convergence.  For  this  model,  convergence  is  determined  by  using  the  sum  of  the 
residuals  squared  as  reference.  This  number  should  approach  zero  as  the  number  of 
elements  is  increased. 

To  test  for  convergence,  executions  of  diatom  were  conducted  for  the  harmonic 
oscillator  in  a  similar  fashion  as  described  earlier  in  the  validation  section.  The  force 
constant  was  set  equal  to  4.0,  and  hJ2n  and  the  reduced  mass  were  set  equal  to  one.  The 
grid  was  constructed  from  -5.0  to  7.0  with  Rc=l  .0  .  Then  diatom  was  executed  to  see 
how  well  it  predicted  the  analytical  eigenvalues  as  measured  by  the  sum  of  residuals 
squared. 

Figure  IV-5  shows  the  results  of  these  executions  as  the  number  of  elements  was 
increased  from  25  to  450.  The  graph  clearly  shows  convergence,  that  is,  as  the  number  of 
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elements  is  increased,  the  solution  approaches  the  analytical  solution  and  the  sum  of  the 
residuals  squared  vanishes. 


Simple  Harmonic  Oscillator:  Sum  of  Residuals 

Squared 


Figure  IV-5.  This  graph  demonstrates  the  convergence  of  the  finite  element  method  for  the  Simple 
Harmonic  Oscillator.  As  the  number  of  elements  is  increased,  the  sum  of  the  residuals  squared 
begins  to  approach  zero.  A  plausible  argument  would  be,  given  an  infinite  number  of  elements  and 
an  infinitely  precise  computer,  the  solution  would  converge  to  the  analytical  solution  exactly.  The 
data  point  at  300  corresponds  to  the  data  presented  earlier  in  the  validation  section. 


Next  the  Morse  Oscillator  was  tested  for  convergence.  Again,  the  executions 
were  set  up  with  De=l 5000.0  and  a=6.0.  The  reduced  mass  and  hJ2%  were  set  to  one. 
The  radius  of  equilibrium  (R^)  was  chosen  as  2.0  .The  left  grid  endpoint  was  set  at  1.883 
and  the  right  endpoint  was  positioned  at  4.500  .  At  this  location,  the  potential  energy 
curve  crosses  the  horizontal  axis.  That  is,  where  the  potential  equaled  zero.  The  thought 
was  that  the  repulsive  branch  above  zero  because  the  solution  of  the  free  states  was  not 
wanted. 
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Figure  IV-6  shows  how  the  solution  did  not  converge  as  expected  when  the 
number  of  elements  was  increased.  The  solution  actually  diverges.  Clearly,  something 
was  wrong  with  the  model.  At  the  time,  the  reasons  were  unclear.  In  fact  at  this  time,  the 
method  of  finite  elements  was  questioned  as  a  reliable  technique  for  these  problems. 


Morse  Oscillator:  Sum  of  Residuals  Squared 
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Figure  IV-6.  This  graph  shows  how  the  solutions  for  the  Morse  potential  did  not  converge.  This 
graph  shows  how  the  solution  diverges  after  the  number  of  elements  is  larger  than  106.  Notice  the 
scale  of  the  Sum  of  Residuals  Squared.  The  residuals  diverge  and  round  off  at  about  104.  This  lack 
of  convergence  is  unacceptable. 


The  code  was  then  tested  thoroughly  because  of  this  divergence.  All  of  the 
algorithms  and  codes  were  combed  through  looking  for  errors.  Then,  the  position  of  the 
right  grid  endpoint  was  tested  to  see  if  it  affected  the  solutions.  For  this  check  the  right 
gridpoint  was  extended  towards  infinity.  This  remedy,  however,  did  not  improve  the 
convergence  problems.  Nothing  seemed  to  work  until  the  left  grid  endpoint  was  moved 
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left  to  1.70  .  This  placement  captured  a  large  part  of  the  repulsive  branch  previously 
thought  unneeded. 

Figure  IV-7  shows  the  solution  for  the  Morse  Oscillator  converging  nicely  to  the 
analytic  solution  with  this  remedy.  The  initial  lack  of  convergence  is  because  the  free 
eigenstates  were  not  being  modeled  correctly. 

The  free  states  were  discovered  to  affect  the  solution  of  the  bound  states.  In 
effect,  when  more  elements  are  added  to  the  solution,  more  eigenstates  are  calculated. 

The  rank  of  the  H  finite  element  matrix  as  discussed  in  Chapter  II  goes  as  2{ne+\),  where 
ne  is  the  number  of  elements.  Since  this  potential  has  only  28  bound  states,  any  ne 
greater  than  13  yields  solutions  beyond  the  last  bound  state.  Therefore,  to  model  these 
free  states  the  positive  repulsive  branch  of  the  potential  must  be  included  in  the 
calculation.  This  requirement  can  be  understood  by  realizing  that  within  the  H  matrix,  a 
lot  of  its  matrix  elements  represent  the  free  states.  When  this  matrix  is  diagonalized,  these 
poorly  represented  states  will  perturb  the  solution  of  the  bound  states. 
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Convergence  of  the  Morse  Oscillator 
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Figure  IV-7.  This  graph  shows  how  the  finite  element  method  does  indeed  converge  for  the  Morse 
Oscillator.  The  left  gridpoint  was  moved  in  from  1.883  to  1.700  capturing  more  of  the  repulsive 
branch.  The  free  states  interact  with  this  part  of  the  potential.  The  solution  of  the  free  states  turns 
out  to  be  important  for  the  convergence  of  the  bound  states. 


Further  testing  explored  the  relationship  between  convergence  of  the  left  grid 
endpoint’s  position.  Figure  IV-8  shows  this  relationship.  As  background,  the  number  of 
elements  for  all  of  these  solutions  was  fixed  at  300,  therefore,  702  eigenstates  were  being 
calculated.  The  solution  can  be  seen  to  slowly  converge  as  the  left  gridpoint  is  moved  to 
the  right.  This  convergence  is  attributed  to  the  shrinking  stepsize.  At  a  critical  point , 
somewhere  around  1.861,  the  solution  diverges.  This  divergence  happens  because  the 
repulsive  branch  of  the  potential  is  truncated  when  the  left  endpoint  is  set  to  the  right  of 
1.861.  The  free  states,  therefore,  won’t  get  modeled  correctly.  These  poorly  formed  free 
states  perturb  the  solution  of  the  bound  states. 

This  discovery  was  the  genesis  of  the  automatic  grid  algorithm.  As  the 
minimization  routine  changes  parameters  of  the  potential,  the  potential’s  repulsive  branch 
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may  swing  off  a  previously  defined  grid.  If  this  situation  happens,  the  solutions  will 
diverge,  and  the  minimization  routine  may  never  find  the  correct  solution.  Therefore,  an 
automatic  system  must  be  in  place  to  ensure  the  potential  is  always  well  defined  on  the 
grid. 


Convergence/Divergence  as  Position  of  Left  Grid  Endpoint 
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Figure  IV-8.  This  graph  shows  the  convergence/divergence  behavior  of  the  finite  element  solution  of 
the  Morse  Oscillator  as  the  left  grid  endpoint  is  moved  about.  The  solution  slowly  converges  as  this 
position  is  moved  from  0  to  1.861 .  This  convergence  can  be  attributed  to  the  decrease  in  the  stepsize. 
The  number  of  elements,  which  was  300,  was  kept  constant  for  of  these  calculations.  The  solution 
then  quickly  diverges  as  the  repulsive  branch  gets  truncated.  The  repulsive  branch  is  required  to 
help  model  the  free  states.  A  poorly  modeled  free  state  perturbs  the  solutions  of  the  bound  states. 


The  solutions  do  converge  given  the  proper  position  of  the  grid  endpoints  have 
been  found.  Assuming  this  is  true,  the  solution  should  get  better  as  the  number  of 
elements  is  increased.  But  at  what  cost?  Recall,  Equation  11-33  established  a  relationship 
between  the  residual  and  the  stepsize.  Taking  the  natural  logarithm  of  this  equation 
yielded 
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ln|  e|xj  j  >  ln(c)  +  4  ln(h) 


(IV-l) 


assuming  cubic  interpolating  polynomials  and  the  residual  is  less  than  one.  The  slope  of 
this  equation  should  be  at  least  4. 

Figure  IV-9  shows  a  convergence  rate  was  far  greater  than  4.  In  fact,  it  averages 
around  10  for  the  stepsizes  tested.  The  data  used  to  generate  this  graph  are  identical  to 
that  used  to  for  Figure  IV-7.  This  curve  was  constructed  by  computing  the  natural 
logarithm  of  each  residual  and  stepsize,  then  the  slope  of  this  curve  was  calculated. 
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Figure  IV-9.  This  graph  shows  the  slope  of  curve  formed  by  evaluating  the  natural  log  of  the 
residual  versus  the  natural  log  of  the  stepsize.  In  this  graph,  the  stepsize  gets  smaller  to  the  right. 
The  slope  of  this  curve  averages  about  10,  far  greater  than  4  as  predicted.  These  calculations  are  for 
the  identical  execution  parameters  as  used  to  generate  Figure  IV-7. 
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However,  there  is  a  cost  for  increasing  the  number  of  elements.  The  rank  of  the 
matrices  to  be  solved  are  2(ne+\),  where  ne  is  the  number  of  elements.  The  cost  for  most 
matrix  solving  algorithms  is  normally  the  rank  of  matrix  cubed. 

To  find  the  actual  cost  of  running  diatom,  the  “user  time”  was  recorded  as  a 
function  of  the  number  of  elements.  The  number  of  elements  was  increased  from  25  to 
450  for  this  measurement, and  the  time  was  measured  on  a  Sun  Sparc20  ®.  Figure  IV- 10 
graphically  shows  these  results  for  both  the  Morse  and  Harmonic  potentials.  The  cost  for 
adding  elements  is  heavy.  In  fact,  anything  over  500  gets  prohibitively  slow. 

This  curve  was  then  found  to  fit  the  polynomial  listed  below 


time  =  8.535  xlO  6  ne3  (seconds)  (IV-2) 

where  time  has  units  of  seconds. 

This  cubic  relationship  correlates  to  the  cost  expected  when  solving  matrices. 

This  cubic  result  implies  that  very  little  time  is  spent  building  the  matrices. 
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ExeaAion  Time  vs  Nuitjer  of  Bemerts 


Number  of  Bements 


Figure  IV-10.  This  graph  illustrates  the  cost  of  adding  elements.  As  seen  previously,  the  more 
elements  used  for  the  solution,  the  better  the  answers  are.  However,  the  time  for  solving  the  finite 
element  method  goes  as  the  number  of  elements  cubed. 
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Hydrogen  Molecule 


In  the  previous  section,  the  finite  element  method  was  shown  to  converge  nicely 
for  the  Harmonic  and  Morse  Oscillators.  But  can  the  full  algorithm,  with  its  non-linear 
minimization  routine,  find  an  unknown  potential? 

To  test  the  diatom  algorithm,  a  numerical  (synthetic)  potential  for  the  H2  (X  *Eg+) 
molecule  was  used  to  calculate  the  vibrational  eigenvalues(l,  2).  This  synthetic  potential 
is  an  accepted  approximation  for  the  H2  potential.  A  variational  technique  was  used  to 
calculate  these  eigenvalues.  This  calculation  included  151  harmonic  wavefunctions  as  a 
finite  basis  set  for  this  approximation  (3).  Due  to  numerical  limitations,  only  the  energy 
values  up  to  state  14  were  used.  Ground  state  H2  has  17  bound  states.  These  eigenvalues 
now  represent  experimental  data,  except  this  synthetic  potential  is  exactly  known  unlike 
that  of  a  molecule. 

The  vibrational  eigenvalues  (in  atomic  units)  were  input  into  diatom  as  the 
experimental  values.  The  reduced  mass  was  input  as  91 1 .422  a.u.,  and  Re  was  set  at 
1 .424  a.u.  The  Morse  potential  was  selected  to  see  if  it  could  model  the  H2  molecule. 

The  program  diatom  had  no  other  information  beyond  this.  A  perfect  fit  was  not 
expected. 

The  program  diatom  was  then  executed  “in  the  blind”  over  six  times.  That  is,  the 
numerical  H2  potential  described  earlier  was  kept  in  confidence  until  all  calculations  were 
completed. 
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After  each  run,  the  parameter  space  was  fine  tuned  as  the  minimization  routine 
began  to  close  in.  At  the  same  time,  the  number  of  elements  was  increased  to  aid  the 
search. 

The  program  diatom  finally  settled  on  an  anharmonicity  term  of  1.0912  inverse 
a.u.  and  a  dissociation  energy  of  0.177  Hartrees.  At  this  point,  the  grid  consisted  of  550 
elements  on  the  interval  0.04  to  13.80  a.u.  An  interesting  note,  these  grid  endpoints  were 
selected  by  the  “smart  grid”  algorithm,  not  by  the  user.  The  reported  dissociation  energy 
is  0.1745  Hartrees  (2:2467),  therefore,  diatom  slightly  over  predicted  this  energy. 

However,  the  new  found  potential  showed  excellent  agreement  with  the  numerical 
H2  potential  energy  curve.  The  Morse  potential,  as  it  turns  out,  provides  for  an  excellent 
approximation  of  the  H2  molecule.  Figure  IV- 1 1  shows  a  Mathematica®  plot  of  these 
two  curves. 
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Figure  IV-11.  This  figure  shows  the  comparison  between  the  potential  diatom  found  for  the  H2 
molecule,  and  the  numerical  potential  used  to  calculate  the  vibrational  eigenvalues.  The  dashed 
curve  is  the  accepted  numerical  H2  molecule  potential,  and  the  solid  curve  is  the  potential  found  by 
diatom.  The  program  diatom  “found  this  curve  in  the  blind.”  The  only  information  it  had  to  work 
with  was  the  reduced  mass  of  the  molecule,  the  equilibrium  position,  and  the  vibrational  eigenvalues. 
The  Morse  function,  therefore,  can  be  used  to  model  the  H2  molecule  with  a  great  deal  of  accuracy, 
though  not  perfect.  The  units  of  this  plot  are  in  atomic  units,  therefore,  the  vertical  axis  is  energy  in 
Hartrees,  and  the  horizontal  axis  is  distance  measured  by  Bohr  radii. 


The  numerical  eigenvalues  also  showed  excellent  agreement.  The  program 
diatom  attempts  to  fit  all  of  the  eigenvalues  simultaneously.  This  process  means  it  tries 
just  as  hard  to  fit  the  upper  states  as  the  lower  states,  but  not  always  with  success.  For  the 
first  1 1  states,  the  relative  error  was  less  than  10%,  and  for  the  bottom  five  levels  it’s  less 
than  2%.  From  this  table,  it  is  clear  the  Hydrogen  molecule  can’t  be  exactly  be  modeled 
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by  the  Morse  function  because  in  some  places  it  under  predicts,  and  others  it  over 
predicts.  The  residuals  do  not  follow  a  systematic  pattern. 


Table  IV-3.  This  table  shows  the  accuracy  of  diatom’s  predictions  for  the  H2  molecule.  The  accepted 
values  are  the  eigenvalues  calculated  using  the  accepted  potential  curve.  The  units  of  these  energies 
are  in  Hartrees. 
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Conclusion 


In  conclusion,  the  program  diatom  and  the  method  of  finite  element  was  proven  to 
provide  credible  solutions  of  Schrodinger's  Equation.  These  solutions  included 
wavefunctions,  energy  eigenvalues,  and  potential  curves  for  diatomic  molecules. 

The  numerical  wavefunctions  of  the  Morse  and  Harmonic  oscillators  were  shown 
to  be  orthonormal.  The  relative  error  of  the  eigenvalues  was  typically  on  the  order  of  10‘6 

•3 

and  10‘  percent  for  the  Morse  and  Harmonic  functions,  respectively.  Convergence  and 
cost  of  convergence  was  also  discussed. 

The  H2  molecule  was  then  used  to  show  the  model  could  be  used  to  find  an 
“unknown”  potential  curve.  The  model  predicted  an  anharmonicity  of  1.0912  inverse  a.u. 
and  a  dissociation  energy  of  0.177  Hartrees  for  the  ground  state  of  H2.  The  graphical 
comparison  of  the  accepted  potential  energy  curve  and  the  predicted  curve  showed 
excellent  agreement. 
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V .  Conclusion 


Summary 

The  method  of  finite  elements  was  proven  to  be  a  viable  technique  for  solving 
Schrodinger's  Equation  for  the  vibrational  states  of  diatomic  molecules.  The  calculation 
of  Franck-Condon  Factors  was  explained  using  experimental  spectroscopic  data  in 
conjunction  with  four  computer  programs  {efit,  dunham,  diatom,  fcf).  The  focus  of  this 
research,  however,  centered  on  the  validation  of  one  of  those  programs,  diatom. 

The  program,  diatom,  requires  as  input  the  experimental  vibrational  eigenvalues 
for  each  electronic  state.  Then  using  the  method  of  finite  elements,  the  program  solves 
Schrodinger's  Equation  and  computes  the  eigenvalues.  These  eigenvalues  are  then 
compared  to  the  experimental  values  to  establish  how  well  the  selected  potential  energy 
model  resembles  the  unknown  real  potential  surface.  The  adequacy  of  this  selected 
potential  model  is  measured  by  the  sum  of  the  residuals  squared.  A  non-linear 
minimization  routine  then  modifies  the  parameters  of  the  potential  energy  model  to 
reduce  this  value,  the  sum  of  the  residuals  squared.  This  process  is  repeated  until  the  sum 
of  residuals  squared  has  been  minimized  identifying  the  optimal  parameter  choice.  After 
the  optimal  potential  model  has  been  isolated,  the  vibrational  eigenfunctions  are  then 
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calculated.  The  vibrational  eigenfunctions  from  two  different  electronic  states  are  then 
used  to  calculate  the  Franck-Condon  Factors. 

To  validate  the  program  diatom,  the  Simple  Harmonic  and  Morse  Anharmonic 
Oscillators  were  exploited  for  their  known  analytical  solutions.  The  relative  error  of  the 
computational  eigenvalues  was  typically  on  the  order  of  1 0  and  1  O'  percent  for  the 
Morse  and  Harmonic  functions,  respectively.  The  numerical  wavefunctions  of  these 
oscillators  were  shown  to  be  orthonormal,  and  matched  the  analytical  wavefunctions 
extremely  well  when  plotted. 

Next,  to  demonstrate  the  model  could  indeed  be  used  to  isolate  an  optimal 
potential  surface,  the  H2  molecule  was  selected  for  the  test.  The  H2  molecule  is  the 
simplest  and  one  of  the  most  studied  molecules,  lending  itself  for  an  this  test  ideally.  The 
vibrational  eigenvalues  for  the  electronic  ground  state  were  generated  using  an  accepted 
potential  surface  for  the  H2  molecule.  These  eigenvalues  were  then  input  into  the  model. 
Using  a  Morse  potential  model,  diatom  isolated  the  optimal  surface  with  an 
anharmonicity  term  of  1 .0912  inverse  a.u.  and  a  dissociation  energy  of  0. 1 77  Hartrees. 
The  graphical  comparison  of  the  accepted  potential  energy  surface  and  predicted  surface 
showed  excellent  agreement. 
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Recommendations 


In  summary,  this  numerical  technique  has  proved  to  be  an  excellent  method  for 
finding  unknown  potential  surfaces  of  real  molecules.  However,  before  this  method  can 
be  shown  to  provide  solutions  which  equal  or  better  the  RKR-IPA  approach,  more  testing 
needs  to  be  accomplished.  The  following  suggestions  outlines  areas  for  additional 
research  and  improvements: 

•  Perform  extensive  testing  of  the  code  on  more  complicated  molecules  to 
demonstrate  the  utility  of  this  numerical  approach.  Preferably,  a  molecule 
such  as  PbO  or  some  other  would  be  chosen  which  has  a  large  empirical 
knowledge  base,  and  has  been  similarly  solved  with  the  RKR  method.  The 
results  of  these  studies  should  be  published. 

•  Validate  the  code  for  a  real  molecule  (such  as  PbO)  at  the  Franck-Condon 
Factor  level.  A  molecule’s  wavefunctions  and  potential  energy  curve  can  not 
be  directly  observed.  Only  the  spectral  locations  and  intensities  of  the 
molecule’s  transitions  can  be  directly  observed.  Experimental  Franck-Condon 
Factors  can  be  obtained  by  removing  the  population  influence  of  the  transition 
intensities,  and  then  normalizing  each  intensity.  The  normalization  factor  is 
the  sum  of  the  intensities  after  the  population  influence  has  been  removed. 
Comparisons  at  this  level  mark  the  code’s  final  test  for  validity. 
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•  Add  more  potential  functions  to  the  menu  the  user  can  choose  from.  This  list 
now  only  includes  the  Simple  Harmonic  Oscillator,  Morse  Anharmonic 
Oscillator,  Lennard-Jones,  and  Mie  potentials.  This  numerical  technique  lends 
itself  very  nicely  to  very  sophisticated  potential  functions  which  may  model 
the  physics  better. 

•  Replace  the  non-linear  minimization  routine  currently  used  with  a  faster  and 
more  efficient  routine.  The  current  routine  is  an  excellent  choice  for  a 
turbulent  parameter  space  with  lots  of  tiny  local  minima.  However,  this 
research  suggest  that  these  parameter  spaces  are  fairly  smooth  and  the 
existence  of  only  one  minimum.  This  enhancement  could  dramatically 
decrease  run  time. 

•  Integrate  all  four  codes  (efit,  dunham,  diatom,  fcf)  into  one  seamless  code. 

This  modification  would  make  the  work  required  of  the  user  simpler.  In  this 
way,  the  user  would  only  have  to  manage  one  input  file.  This  enhancement  is 
necessary  before  this  code  is  distributed  in  any  way. 

•  Consider  replacing  the  IMSL®  math  routines  with  non-proprietary  code.  This 
modification  would  eliminate  any  portability  issues  which  may  arise. 
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